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THEORETICAL  STUDIES  AND  DATA  ANALYSIS  OF 
WAVE  PROPAGATION  IN  RANDOM  MEDIA 


FINAL  TECHNICAL  REPORT 
December  1,  1983  to  September  30,  1985 
Contract  No.  MDA903-83-C-0515 
Principal  Investigator:  Stantey  M.  Flatte' 

SUMMARY 

\ 

The  objective  of  this  work  was  to  develop  a  general  theoretical  framework  for  cal¬ 
culating  fluctuations  of  signals  on  waves  propagated  through  random  media  (WPRM) 
and  to  apply  this  framework  to  sound  through  the  ocean;  light  through  the  atmosphere; 
radio  waves  through  the  ionosphere,  solar  wind,  or  interstellar  plasma;  and  any  other 
similar  case  of  waves  propagating  through  continuous  media.  Comparison  with  real 
data  is  an  important  aspect  of  the  effort. 

This  report  will  consist  of  a  summary  {with  list  of  references),  followed  by  copies  of 
the  journal  articles  resulting  from  this  contract  that  have  been  published,  or  have  been 
submitted  for  publication.  Some  work  in  progress  will  not  be  included  in  detail  here, 
since  it  will  be  discussed  in  the  reports  of  future  DARPA  contracts. 

The  two  most  common  signals  sent  on  a  carrier  are  the  phase  and  amplitude  of  a 
nearly  monochromatic  wave.  If  enough  bandwidth  is  available,  one  can  send  a  pulse, 
and  one  can  speak  of  the  intensity  and  arrival  time  of  that  pulse.  The  technical  prob¬ 
lem  is  then  to  explain  the  statistical  behavior  of  the  intensity  and  arrival  time  in  terms 
of  medium  fluctuations,  where  the  medium  is  described  statistically,  usually  by  a  power 
spectrum  covering  a  large  dynamic  range  of  scales.  ~<2 

4 

The  eventual  practical  applications  of  an  understanding  of  WPRM  to  science  and 
to  the  defense  department  are  myriad.  The  phase  of  a  light  wave  from  an  astronomical 
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object  or  a  satellite  is  used  by  a  telescope  to  focus  to  a  detector;  the  quality  of  the 
focus  depends  on  the  state  of  the  atmosphere  •  a  random  medium.  A  ground-based 
laser  with  large  optics  attempts  to  focus  on  a  small  spot  for  a  period  of  time;  the  atmo¬ 
sphere  spreads  out  the  spot  by  its  action  on  the  phase  of  the  wave.  Determination  of 
pulsar  parameters  depends  on  observations  of  radio  pulses  through  a  distorting  random 
medium-interstellar  plasma.  Communication  with  spacecraft  by  radio  pulses  using 
radio  telescopes  depends  on  coding  schemes  and  antenna  control  that  must  contend 
with  effects  due  to  the  solar  wind  and  the  earth’s  ionosphere.  Communication  with 
earth  satellites  must  contend  with  the  ionosphere,  and  the  effects  of  a  disturbed  iono¬ 
sphere  (due,  e.g.,  to  nuclear  explosions)  must  be  predicted.  Probing  ocean  processes, 
from  large-scales  (Gulf  stream)  on  down  (internal  waves  and  microstructure)  depends  on 
understanding  acoustic  propagation  through  a  random  medium;  for  example,  a  favored 
method  is  to  send  pulses  over  long  range  and  observe  their  arrival-time  variations. 
Detection  of  submarines  by  passive  acoustics  is  limited  by  ocean  fluctuations  that  con¬ 
trol  the  maximum  antenna  size  and  integration  time  that  can  be  used  for  coherent  sig¬ 
nal  integration.  Active  sonar  for  communication  has  a  limited  bandwidth  due  to  ocean 
fluctuations.  Determination  of  the  characteristics  of  earthquakes,  thought  to  be  impor¬ 
tant  for  earthquake  prediction,  is  done  solely  by  seismic  detection.  The  earth  through 
which  the  seismic  waves  travel  has  randomness,  and  this  limits  the  information  that 
can  be  gleaned  from  seismic  signals.  The  same  effect  limits  our  ability  to  distinguish 
between  underground  nuclear  explosions  and  earthquakes,  and  thus  affects  our  political 
stance  vis-a-vis  nuclear  test  bans.  On  the  other  band,  the  distortions  of  seismic  waves 
due  to  the  earth’s  random  properties  can  be  used  as  a  probe  of  those  properties  and 
hence  can  lead  to  a  better  understanding  of  earth  structure. 

The  path-integral  method  for  treating  wave  propagation  has  been  successfully 
used  by  the  principal  investigator  for  the  analysis  of  many  experiments  in  ocean  acous¬ 
tics.1  This  method  is  therefore  utilized  extensively  in  our  work.  Another  school  of  wave 


propagation  theory,  which  began  in  the  early  1960*s  in  order  to  explain  experiments  in 
light  transmission  through  the  turbulent  atmosphere,  utilizes  a  different  technique 

A 

involving  partial  differential  equations  for  the  moments  of  the  wave  field.  The  moment 
equations  and  the  path  integral  have  now  been  used  by  enough  researchers  that  the 
value  of  both  approaches  is  appreciated,  but  the  relations  between  the  two  methods  has 
remained  confusing  to  many  people. 

A  first  step  in  WPRM  is  characterizing  the  statistics  of  the  random  medium.  In 
this  report  we  will  speak  of  the  spectrum  of  medium  fluctuations,  which  would  be 
obtained  by  dragging  an  index-of-refraction  sensor  through  the  medium  and  taking  a 
Fourier  transform  of  the  resulting  time  (=  space)  series.  This  spectrum  is  characterized 
by  a  power  law  (e.g.  -5/3  for  Kolmogorov  turbulence)  that  implies  much  more  variation 
at  large  scales  than  at  small  scales.  It  is  likely  that  a  sensor  dragged  in  different  direc¬ 
tions  will  observe  different  spectra  even  on  the  average.  In  that  case  we  speak  of  an 
anisotropic  spectrum.  For  example,  the  strength  of  the  ocean  spectrum  is  much  higher 
in  the  vertical  than  the  horizontal  for  the  same  wave-number.  Finally  it  is  important 
to  know  that  the  spectrum  cannot  continue  indefinitely  at  either  large  or  small  scales. 
At  the  "outer  scale"  the  spectrum  cuts  off  due  usually  to  finite  container  size  —  the 
height  of  the  atmosphere  or  the  depth  of  the  ocean.  At  the  "inner  scale"  the  spectrum 
cuts  off  due  to  physical  processes;  for  example,  viscosity  becomes  important  at  scales  of 
order  a  few  millimeters  in  the  atmosphere  and  ocean. 

This  technical  report  covers  the  two  year  period  of  our  contract  effort.  The  next 
several  paragraphs  summarize  the  technical  results  we  have  obtained:  more  details  are 
given  in  following  sections.  The  work  has  been  carried  out  under  the  direction  of  Dr. 
Stanley  Flatte',  and  involves  effort  by  Dr.  Flatte',  senior  scientist  Dr.  Frank  Henyey, 
two  post-doctoral  researchers  Drs.  Dennis  Creamer  and  Rod  Freblicb,  and  a  graduate 
student  (in  the  UCSD  Electrical  Engineering  and  Computer  Science  Department) 
Johanan  Codona. 


Our  progress  in  understanding  the  travel  time  or  pulses  in  random  media  has 
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resulted  in  a  paper  published  in  Physical  Review  Letters,  as  well  as  results  that  will 
lead  to  later  publication.  Pulses  sent  through  a  fluctuating  medium  arrive  earlier  or 
later  than  they  would  in  the  absence  of  fluctuations,  depending  on  the  particular  reali¬ 
zation  of  the  medium.  The  variance  of  arrival  time  can  be  calculated  by  straightfor¬ 
ward  methods  in  the  geometrical  optics  limit.  Our  dramatic  new  result  is  for  the  aver¬ 
age  arrival  time,  which  we  find  advanced  in  weak  media.  Heretofore,  researchers  were 
of  the  opinion  that  pulses  were  delayed  on  the  average.  This  effect  is  of  little  impor¬ 
tance  in  communication  applications  where  the  average  arrival  time  is  usually  less 
important  than  the  variance.  However,  the  effect  can  be  important  in  probing  a  ran¬ 
dom  medium  for  large-scale  variations  by  their  effect  on  average  travel  time.  Our 
result  implies  a  possible  confusion  between  a  changing  turbulence  level  and  a  change  in 
the  average  index  of  refraction  on  a  large  scale.  For  example,  ocean  acoustic  tomogra¬ 
phy  attempts  to  measure  the  warming  of  a  100-km-square  area  of  the  ocean  by  an 
expected  change  in  travel  time  of  about  20  ms.  However  we  find  a  change  in  average 
travel  time  of  about  10  ms,  due  to  an  internal-wavefield  that  has  no  average  warming 
at  all.  We  have  recently  studied  the  range  dependence  of  this  effect,  and  have  found 
that  it  grows  as  the  square  of  the  range.  This  implies  that  experiments  being  planned 
in  the  1000-4000  km  region  will  have  major  difficulties  sorting  out  the  effects  of  internal 
waves  from  the  effects  of  large-scale  structure.  Most  importantly,  the  determinations  of 
internal-wave  effects  will  NOT  be  contaminated  by  the  large-scale  effects.4 

The  above  understanding  of  travel-time  effects  arose  from  studying  the  mutual 
coherence  function  (MCF)  of  the  complex  wave  function  of  the  wave  field  arriving  at  a 
receiver.  We  have  developed  quantitative  treatments  of  the  MCF  in  an  anisotropic 
medium  with  curved  deterministic  rays,5  and  have  applied  these  treatments  to  data 

c 

from  a  35-km,  5-kHz  ocean-acoustic  experiment,  with  good  success.  The  medium 
fluctuations  in  that  experiment  were  measured  by  instruments  that  were  independent  of 
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the  acoustic  information,  so  no  free  parameters  were  available  to  the  theory. 

Moving  from  arrival  time,  which  is  related  to  phase,  we  discuss  amplitude  or  inten¬ 
sity.  We  have  developed  a  method  for  calculating  the  spatial  correlation  function  of 
intensity  on  a  transverse  plane  through  a  receiver.  This  is  a  long-standing  problem 
that  is  of  great  importance  in  using  wave  propagation  for  probing  the  structure  of  a 
random  medium,  because  measuring  intensity  is  often  the  observation  that  can  be  made 
most  easily.  In  addition,  an  amplitude-modulated  signal  will  be  degraded  by  intensity 
fluctuations  due  to  the  medium.  The  standard  theory  develops  a  series  solution  for  the 
intensity  spatial  spectrum.  The  first  few  terms  are  an  accurate  representation  of  the 
small-wave-number  end  of  the  spectrum.  In  order  to  calculate  the  high-wave-number 
region  many  terms  of  the  series  had  to  be  evaluated.  We  have  determined  a  different 
series  expansion,  whose  first  few  terms  give  the  high-wave-number  section  of  the  spec¬ 
trum.  Hence  the  evaluation  of  the  full  spectrum  is  simplified  considerably.  We  have 
also  made  considerable  progress  toward  evaluating  the  intensity  spectrum  for  an  arbi¬ 
trary  source  distribution,  going  beyond  the  standard  procedure  of  considering  the  spe¬ 
cial  cases  of  a  point  source  or  an  incident  plane  wave.  Our  general  case  will  include  a 
source  that  is  extended  over  a  large  aperture.  An  example  of  a  coherent  source  of  large 
aperture  would  be  a  large-aperture  laser  beam.  An  example  of  an  incoherent  source  is  a 

planet,  or  an  illuminated  satellite,  or  an  infrared  plume  from  an  ascending  booster.  A 

7 

paper  describing  our  results  is  in  review  at  the  journal  of  Radio  Science. 

Intensity  correlations  at  two  different  frequencies  are  of  interest  for  a  variety  of 
reasons.  We  have  derived  the  intensity  cross-spectrum  for  scintillations  caused  by  a 
plane  wave  passing  through  a  random  phase  screen.  A  common  approximation  for  a 
case  of  this  sort  is  the  Gaussian-Field  approximation,  in  which  the  cross-spectrum  is 
modelled  as  the  transform  of  the  square  of  the  second  moment.  We  have  shown  how 
this  approximation  breaks  down  when  the  outer  scale  is  large  compared  with  the  diame¬ 
ter  of  the  scattering  disk  (the  transverse  region  of  significant  wave  energy).  A  paper 
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describing  these  results  has  been  submitted  to  Radio  Science.  Furthermore,  the  thesis 
of  J.L.  Codon  a  will  contain  much  of  the  two-frequency,  extended-source  results.9 

Because  researchers  favoring  the  moment-equation  method  or  the  path-integral 
method  typically  knew  only  one  of  the  methods  in  any  depth,  the  relation  between  the 
two  methods  has  been  a  mystery  to  many.  We  have  expended  quite  a  bit  of  effort  to 
understand  this  relation.  We  have  shown  that  the  two  methods  are  mathematically 
equivalent,  in  much  the  same  way  that  the  Heisenberg  and  Schrodinger  approaches  to 
quantum  mechanics  were  shown  to  be  mathematically  equivalent.  A  better  analogy  for 
those  familiar  with  quantum  mechanics  is  the  equivalence  of  the  Schrodinger  and  Feyn¬ 
man  approaches  to  quantum  mechanics.  The  equivalence  extends  to  the  equivalence 
term  by  term  of  the  series  solutions  for  the  intensity  spectra  mentioned  earlier,  a  paper 
describing  these  results  has  been  accepted  for  publication  in  the  Journal  of  Mathemati¬ 
cal  Physics.10 

In  nearly  all  cases,  in  order  to  compare  theory  to  experimental  data  in  WPRM,  we 
must  use  a  model  spectrum  for  the  medium  fluctuations.  We  have  developed 
phenomenological  spectra,  as  a  function  of  wave  vector,  that  allow  for  an  anisotropic 
component  added  to  a  turbulent  isotropic  component.11  This  model  is  meaningful  both 
for  the  ocean,  where  the  anisotropic  component  represents  internal  waves,  and  the  iono¬ 
sphere,  where  the  anisotropy  is  due  to  electrons  preferentially  moving  along  magnetic 
field  lines.  We  are  in  the  process  of  calculating  intensity  spectra  in  the  weak  fluctua¬ 
tion  region  using  these  model  spectra.  We  have  data  from  an  ocean-acoustic  experi¬ 
ment  that  will  be  used  for  comparison  purposes;  the  experiment  utilized  10-70  kHz 
sound  over  several  hundred  meters  under  the  Arctic  ice.  We  should  note  that  the 
weak-fluctuation  regime  is  one  in  which  the  intensity  series  solution  for  the  low-wave¬ 
number  regime  is  the  only  relevant  one. 

For  more  than  one  hundred  years,  eclipse  observers  have  noted  "mysterious" 
bands  of  shadows  moving  on  the  ground  just  before  and  after  an  eclipse.  Many  exotic 
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theories  of  these  shadow  bands  have  been  put  forward,  but  most  observers  agree  that 
they  are  probably  due  to  atmospheric  scintillation  that  becomes  visible  when  the  cres¬ 
cent  of  the  eclipsed  moon  becomes  thin  enough  (a  few  minutes  before  and  after  an 
eclipse).  Johanan  Codona,  the  graduate  student  associated  with  this  project,  has  made 

i  o 

the  first  systematic  application  of  WPRM  theory  to  eclipse  shadow-band  observations. 
He  explains  the  orientation  and  contrast  of  the  bands  as  a  function  of  time,  and 
describes  the  effects  of  eclipse  geometry  and  the  importance  of  wind  direction.  He 
relates  shadow-band  observations  to  stellar-scintillation  observations.  One  important 
conclusion  he  draws  is  that  as  the  illuminated  crescent  gets  thinner,  the  shadow-band 
observations  probe  higher  into  the  atmosphere.  Recently  published  data14  from  the 
eclipse  of  February  16,  1980  in  India  agrees  with  Codona’s  predictions.  Further  data 
from  the  annular  eclipse  of  May  30,  1984  in  Georgia  should  soon  be  forthcoming. 

We  have  begun  the  analysis  of  seismic  data  from  the  Center  for  Seismic  Studies. 
Two  nuclear  explosions  in  the  Soviet  Union  with  good  detections  on  the  NORESS  array, 
which  has  about  twenty  elements  spaced  out  to  a  few  kilometers  have  been  obtained. 
The  first  look  shows  rather  small  travel-time  fluctuations,  somewhat  at  odds  with  the 
large  amplitude  fluctuations  that  have  been  suggested  previously.  Some  of  the  prelim¬ 
inary  data  are  shown  in  a  subsequent  section.  We  have  looked  at  the  data  as  a  func¬ 
tion  of  frequency  up  to  about  20  Hz,  and  have  seen  no  obvious  systematic  differences  in 
travel  time  between  the  different  frequencies,  except  for  an  unusual  change  in  the 
arrival  structure  between  5  and  10  Hz. 

The  theory  of  wave  propagation  through  three-dimensional,  continuous,  random 
media  is  based  largely  on  the  parabolic  wave  equation.  This  equation  did  not  appear  in 
a  classical  physics  context  until  about  1950.  Yet,  after  1926  it  was  used  in  a  quantum- 
mechanical  context,  where  it  is  called  the  Schrodinger  Equation.  We  have  studied  the 
history  of  the  development  of  this  equation,  and  have  developed  some  pedagogically  use¬ 
ful  ideas  on  how  to  introduce  the  Schrodinger  equation  to  students  in  a  manner  that  is 


less  mysterious  than  is  used  in  the  present  curriculum.  A  paper  on  this  subject  has  been 
submitted  to  the  American  Journal  of  Physics.15 

We  have  presented  selected  portions  of  our  work  at  meetings  of  the  American  Phy¬ 
sical  Society  at  Providence  (November,  1984),  the  Acoustical  Society  of  America  at  Aus¬ 
tin  (April,  1985),  and  the  Union  of  Radio  Scientists  at  Vancouver  (June  1985),  as  well  as 
in  an  invited  talk  at  the  International  Symposium/Workshop  on  Multiple  Scattering  of 
Waves  in  Random  Media  and  Random  Rough  Surfaces  at  the  Pennsylvania  State 
University  (July,  1985). 16 

Our  longer-term  goals  will  include  the  implementation  of  our  new  theoretical 
results  into  computer  codes  for  calculation  of  general  phase  and  amplitude  fluctuations. 
Two  directions  are  contemplated  that  will  require  large-scale  computing.  The  first  is 
propagating  waves  through  individual  realizations  of  random  media  to  compare  with 
our  theoretical  results  and  to  extrapolate  those  results  into  parameter  regimes  in  which 
the  theory  is  not  valid.  This  propagation  can  be  done  via  a  parabolic  wave  equation,  so 
that  it  is  a  marching  solution.  The  second  involves  evaluating  the  theoretical  formulas 
which  involve  either  multidimensional  ordinary  integrals,  or  in  some  cases,  path 
integrals.  We  have  begun  some  computer  work  on  simulation  of  WPRM  using  our  VAX, 
and  we  plan  to  implement  the  code  on  an  IBM-PC  that  has  two  FFT  hardware  boards 
that  should  allow  uninterrupted  calculations  at  about  twice  the  speed  of  a  VAX.  These 
efforts  are  in  preparation  for  proposals  to  do  simulation  calculations  on  a  CRAY. 

Finally,  it  is  desirable  to  make  the  comparison  in  a  unified  way  between  these 
theoretical  approaches  and  data  from  experiments  in  seismology,  ocean  acoustics, 
atmospheric  optics,  and  radio  waves  through  the  ionosphere,  the  solar  wind,  or  the 
interstellar  medium.  Our  results  will  be  important  building  blocks  in  making  that  hope 
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It  is  pointed  out  that  a  continuous  random  medium  can  cause  an  average  advance  of  the  arrival 
time  of  a  pulse.  This  advance  will  occur  for  unsaturated  and  partially  saturated  propagation,  but  not 
in  full  saturation  (which  corresponds  to  the  discreie-scatterer  case).  The  effect,  which  is  associated 
with  Fermat’s  principle  of  least  time,  can  be  observed  by  measuring  the  difference  between 
intensity-weighted  and  unweighted  average  arrival  times. 

PACS  numbers:  03.40.Kf,  05.40.  +j.  42  20.Cc,  43  30  +  m 


Wave  packets,  or  pulses,  are  frequently  used  to 
probe  inhomogeneous  media.  If  the  wave  speed  in  the 
medium  varies  on  scales  small  compared  with  the  total 
distance  traveled  by  the  pulse,  then  the  medium  is 
treated  by  statistical  methods.  “Macroscopic"  exam¬ 
ples  include  sound  through  ocean  internal  waves,1,3 
light  through  atmospheric  turbulence,3  and  radio 
waves  through  plasmas  such  as  the  ionosphere,4  the 
solar  wind,5  or  the  interstellar  medium.6  Microscopic 
examples  include  various  “sounds”  through  liquid 
helium,7  and  waves  through  inhomogeneous  con¬ 
densed  matter.*  These  continuous  media  may  be  dis¬ 
tinguished  from  media  consisting  of  discrete  scatterers 
such  as  would  occur  in  light  transmission  through  fog,9 
or  in  wave  transmission  through  a  gaslike  medium 
with  random-point  particles.10 

This  Letter  points  out  that  a  fluctuating  continuous 
medium  can  cause  an  average  advance  of  the  pulse  ar¬ 
rival  time.  All  previous  analyses  have  dealt  with  situa¬ 
tions  in  which  puises  are  delayed  on  the  average.1113 
By  convention,1*13  the  ensemble  average  of  a  random 
medium  is  taken  as  the  medium  reference  state,  and 
the  small  fluctuations  about  this  reference  state  are 
thus  by  definition  a  zero-mean  random  process.  The 
arrival-time  advance  or  delay  is  relative  to  the  travel 
time  through  the  reference  states.  Thus,  for  example, 
results  through  turbulent  air  or  plasma  are  relative  to 
quiescent  air  or  plasma,  not  vacuum. 

The  behavior  of  a  wave  propagating  through  a  ran¬ 
dom  medium  is  controlled  by  relationships  between 
the  wave  number  <*)  of  (he  propagating  wave,  the 
range  (/?),  and  the  strength  and  size  of  the  medium 
fluctuations.1,3  Unsaturaied  behavior  corresponds  to 
one  stationary-phase  path  (ray),  and  occurs  if  the 
medium  fluctuations  are  weak  enough.  In  fully  saturat¬ 
ed  behavior  the  original  ray  breaks  up  into  many  new 
microrays  which  are  statistically  independent  of  each 
other.  Propagation  through  a  medium  of  discrete 
scatterers  falls  in  this  category.  Partially  saturated 
behavior  occurs  in  a  strongly  fluctuating  medium  with 
a  power-law  spectrum,  which  has  enough  small-scale 


fluctuations  to  cause  the  breakup  into  many  microrays, 
and  enough  large-scale  fluctuations  to  make  the  mi¬ 
croray  bundle  behave  like  a  single  ray  in  its  wandering 
from  the  unperturbed  ray.  Experiments  in  waves 
propagating  through  continuous  random  media  typical¬ 
ly  fall  into  this  category.  We  deal  only  with  the  impor¬ 
tant  case  in  which  the  transverse  wandering  from  the 
unperturbed  ray  is  small  compared  with  the  range  of 
propagation. 

Briefly  our  results  are  as  follows:  If  the  travel  time 
of  a  pulse  is  averaged  over  an  ensemble  of  the  random 
medium,  with  each  pulse  weighted  by  its  intensity, 
then  the  average  pulse  is  delayed,  regardless  of  the 
type  of  propagation  behavior,  in  agreement  with  previ¬ 
ous  results.11,13  However,  if  the  average  travel  time  is 
obtained  without  weighting  by  pulse  intensity,  then  a 
pulse  advance  is  expected  for  both  unsaturated  and 
partially  saturated  behavior,  while  a  pulse  delay 
remains  for  the  fully  saturated  case.  The  difference 
between  intensity-weighted  and  unweighted  travel 
time  probes  the  variance  of  the  first  derivative  of  the 
refractive  index,  smoothed  over  a  microray  bundle. 

To  explain  our  effect  qualitatively  we  first  take  a 
simple  special  case.  Consider  a  point  source  and  point 
receiver  separated  by  range  R,  and  a  homogeneous 
medium  in  the  absence  of  fluctuations,  so  that  the  un¬ 
perturbed  ray  from  source  to  receiver  is  a  straight  line. 
The  random  medium  is  concentrated  in  a  “phase 
screen"  at  a  distance  z  from  the  source.  This  screen 
has  the  effect  of  advancing  the  time  of  a  wave  front  by 
a  random  amount  tlx)  where  oris  the  position  on  the 
screen,  and  /(x)  is  a  stationary  Gaussian  random  pro¬ 
cess  with  zero  mean.  (We  take  x  as  one-dimensional 
for  simplicity.) 

Weak  fluctuations.— In  the  geometrical-optics  limit 
only  one  ray  exists  from  source  to  receiver.  The  travel 
time  for  a  path  through  point  x  is 

Tix)  =*  T0+O.Sco]Ax2- r(x),  (1) 

where  A~'mz(R-z)/R.  By  Fermat’s  principle  the 
ray  is  at  the  point  x,  such  that  Tlx,)  is  a  minimum. 


©  1985  The  American  Physical  Society 


9 


Volume  55,  Number  ! 


PHYSICAL  REVIEW  LETTERS 


I  July  1915 


For  the  case  of  weak  fluctuations  we  may  expand  Hx) 

as 

i(x)-/o+^  +  0.5rxJ.  (2) 

The  position  of  the  ray  follows  to  first  order  as 

x,-i4"'cor'.  (3) 

The  travel  time  of  the  ray  is  then 

T(x,)  -  70+  0.5<v< "  V*-  r0-  c0A  '  V1.  (4) 

This  case  requires  that,  typically 

(5) 

But  r  (and  hence  r0  and  /')  are  (by  construction)  ran* 
dom  variables  with  zero  mean.  Therefore  the  /<>  term 
will  disappear  in  the  average  travel  time  and  the  only 
effect  of  the  fluctuations  will  come  from  the  t‘  terms. 
These  terms  arise  because  the  ray  has  moved  away 
from  its  unperturbed  position.  The  first  /'  term  is  pos¬ 
itive,  corresponding  to  a  pulse  delay,  and  represents 
the  effects  of  geometry;  the  perturbed  path  is  physical* 
ly  longer  than  the  unperturbed  one.  The  second  r‘ 
term  is  negative,  corresponding  to  a  pulse  advance;  we 
cal)  this  the  Fermat  term;  the  ray  sought  out  a  region 
of  the  medium  with  a  pulse  advance.  The  Fermat 
term  is  twice  as  large  in  magnitude  as  the  geometry 
term.  The  average  travel  time  is 

W-To-OScoA-'U*).  (6) 

so  that  the  pulse  on  the  average  arrives  early. 

There  is  a  subtlety  to  this  result.  In  the  weak* 
fluctuation  limit  the  intensity  is  controlled  by  the 
focusing  due  to  the  curvature  of  the  wave  front  as  it 
exists  from  the  phase  screen.  It  is  not  difficult  to  show 
that  the  intensity  /  is,  to  first  order, 

/-l+v4_,c0r".  (7) 

Consider  the  intensity-weighted  average  travel  time 

(/ru)>-r,-o.5c(^-,</'J> 

-coA-'U oO.  <8) 

where  the  last  term  comes  from  the  correlation 
between  the  intensity  and  the  travel  time.  For  any 
random  function  /(x)  whose  Fourier  components  are 


uncorrelated  (i.e.,  the  correlation 
translation-invariant) 

function  is 

</eO --(f*). 

(9) 

Therefore 

(mx,)>-r#+0.5r(>y4-'(r,,>. 

(10) 

In  other  words,  the  intensity-weighted  average  travel 
time  is  delayed  by  fluctuations  by  exactly  the  amount 
that  the  unweighted  average  is  advanced!  The  focus¬ 


ing  effect  exactly  cancels  the  Fermat  term,  leaving  a 
resultant  equal  to  the  geometry  effect  alone.  This  oc¬ 
curs  because  a  positive  fluctuation,  which  delays  the 
pulse,  acts  as  a  converging  lens  to  increase  the  inten¬ 
sity. 

The  simple  example  of  a  phase  screen  in  the  weak- 
fluctuation  geometrical-optics  limit  has  illustrated  our 
point.  We  will  now  make  some  remarks  on  generaliza¬ 
tions  to  extended  media  and  strong  fluctuations  which 
we  have  treated  rigorously  but  do  not  have  space 
within  the  Letter  format  to  describe  in  detail.  We  then 
describe  a  rigorous  extension  of  these  results  by  means 
of  a  path-integral  method  to  include  diffractive  effects 
in  a  power-law  medium. 

There  is  no  difficulty  in  extending  the  above  results 
from  a  phase  screen  to  extended  media  in  which  (6) 
and  (10)  are  replaced  by 

<r)-r0--0.5co“,/dr^-,(z) 

xl/*'p*(z.z)l.  (11) 

(/r>-r0-  +  o.sco-’/dz  a  -Hz) 

xl/«fcWz.z')J.  (12) 

p„(z,z')“  (3,m(z)0xjz (*')).  (13) 

where  dx/i(z)  is  the  transverse  gradient  of  the  refrac¬ 
tive  index  due  to  the  fluctuations  at  location  z  along 
the  unperturbed  ray.  These  results  require  the  Mar¬ 
kov  approximation  [that  is,  the  quantity  in  square 
brackets  in  (11)  is  a  local  function  of  zl.  If  an  incident 
plane  wave  rather  than  a  point  source  is  used,  all  three 
terms  (geometry,  Fermat,  and  focusing)  are  reduced 
by  a  factor  of  3.  If  the  Markov  approximation  is  not 
made,  the  ratio  between  Fermat  and  geometry  remains 
-  2,  while  all  terms  are  modified  by  terms  of  order 
Lf/R,  where  Lf  is  the  longitudinal  correlation  length 
of  the  medium  fluctuations. 

Strong  fluctuations.— If  the  medium  fluctuations  are 
strong  enough,  one  can  show  in  the  limit  of  small 
wavelength  that  both  the  Fermat  and  focusing  terms 
become  negligible.  This  occurs  when  the  intensity 
fluctuations  become  of  order  unity.  From  (7)  this  cor¬ 
responds  to 

A-'coir^M^l.  (14) 

If  this  condition  is  satisfied,  the  unperturbed  ray 
breaks  up  into  many  microrays  that  are  not  minima, 
but  extrema,  in  accord  with  Hamilton's  principle  of 
stationary  action.  The  converging  (or  diverging) 
lenses  now  are  so  strong  that  caustics  occur,  destroying 
the  correlation  between  intensity  and  pulse  delay.  The 
Fermat  term  disappears  because  the  microrays  are  not 
at  global  minima;  hence  the  geometry  term  strongly 
dominates,  resulting  in  the  validity  of  (10)  even  for 
unweighted  average  travel  time. 
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If  t(x)  has  a  power-law  spectrum,  then  the  above 
argument  must  be  modified  to  separate  the  effects  of 
large-scale  and  small-scale  fluctuations.  The  small- 
scale  fluctuations,  if  they  alone  satisfy  (1 4).  break  the 
unperturbed  ray  into  a  bundle  of  microrays  of  size  L 
Fluctuations  larger  than  correlate  the  travel  times 
of  all  the  microrays  as  though  they  were  a  single  ray. 
and  cause  the  average  pulse  arrival  time  to  behave  ac¬ 
cording  to  the  weak-fluctuation  rule  (6).  where  </'J)  is 
to  be  interpreted  as  coming  only  from  scales  larger 
than  Lp.  It  then  becomes  crucial  to  estimate  For 
the  phase-screen  case,  we  find 

L„  =  <V4-V'I>,/2.  (15) 

The  difficulty  here  is  defining  the  average  travel  time 
of  a  pulse  that  itself  has  complicated  structure  due  to 
microrays. 

Path-integral  result.— A  major  limitation  of  the 
above  treatment  is  its  restriction  to  cases  of  very  small 
wavelength,  and  the  related  requirement  of  defining 
an  inner  scale  for  the  medium  fluctuations  (in  order 
for  (r"2)  or  even  (/'2)  to  be  finite).  The  path-integral 
method  can  treat  the  case  of  finite  wavelength  for 
intensity-weighted  travel  time. 

To  arrive  at  average  pulse  travel  time  we  note  that 

{IT) -T0 

*=  -  (i/co)9«(|J'*(k0+  k  )>]j(k0))  |*«0.  (16) 

where  k 0  is  a  central  wave  number  of  the  propagating 
wave,  k  is  a  deviation  wave  number  whose  excursion 
represents  the  pulse  bandwidth,  and  0  is  the  reduced 
wave  function  at  the  receiver.  Let  us  treat  the  phase- 
screen  case  for  simplicity.  The  second  moment  as  a 
function  of  wave  number  is  known  to  be’ 3,14 

<«|»’(*0+*).J»Uo>> 

**expl  —  0.5 (f$))C(A),  (17) 

Q{k)~  N  j du  exp\\\iAk$k~'u1  -  k$c$D(u)]\. 

(18) 

where  D(u)  is  the  structure  function  of  the  medium, 

D{u) «  ((/(«/)  -  / (0)  12> .  09) 

and  N  is  a  normalization  such  that  Q( 0)  —  1.  It  is  easy 
to  see  that  Q(k)  has  an  increasing  imaginary  part  as  k 
grows  from  zero. 

At  very  small  k,  Q(k )  is  controlled  by  the  behavior 
of  D(u)  at  small  u,  which  must  be  quadratic  if  an 
inner  scale  exists.  Thus  for  small  u 

D{u)-{t'*)u'.  (20) 

Then’4 

<?(*)-  11-W'C,</'2>*1",/2.  (2D 


and  the  result  for  the  intensity-weighted  average  travel 
time  is  (10)  exactly.15 

If  we  attempt  to  model  D(u)  at  infinitesimal  i/as  a 
fractional  power  law.  then  the  average  travel  time 
diverges. 

This  result  for  finite  wave  number  predicts  a  delay 
equal  to  the  geometry  term.  Most  importantly  this 
derivation  has  made  no  distinction  between  the  weak- 
and  strong-fluctuation  cases,  and  can  be  easily  general¬ 
ized  to  the  extended  medium.  We  can  conclude  that 
{IT)  -  To  is  equal  to  the  geometry  term  alone  regard¬ 
less  of  the  fluctuation  strength,  a  result  already  sug¬ 
gested  by  the  geometrical-optics  calculation. 

When  (18)  is  generalized  to  an  extended  medium  by 
a  path-integral  formalism  the  same  conclusions  are 
easily  shown,  subject  to  the  additional  assumption  of 
the  Markov  approximation. 

An  important  modification  of  the  above  result  oc¬ 
curs  if,  in  the  absence  of  fluctuations,  the  medium  has 
focusing  properties.  In  ocean  sound  propagation  this  is 
due  to  the  sound  channel.  In  radio- wave  propagation 
from  pulsars  this  might  be  due  to  very  large-scale 
medium  fluctuations  that  are  effectively  frozen  during 
the  time  of  observation.  The  modification  can  be  sim¬ 
ply  expressed  by  generalizing  A  "’(z)  in  (11)  and  (12) 
for  a  curved  ray.2  The  key  result  is  that  A  ~  ’(z)  can  be 
negative  for  various  regions  along  the  ray,  and  hence 
the  geometry  term  can  be  negative  for  curved  rays.  This 
complication  is  crucial  to  the  comparison  between  cal¬ 
culation  and  experiment  in  the  ocean,  though  probably 
not  in  other  media.  Note  that  this  effect  provides 
another,  different  mechanism  by  which  fluctuations  in 
a  medium  may  cause  an  average  pulse  advance.  Final¬ 
ly,  dispersive  propagation,  such  as  occurs  for  radio 
waves  through  plasma,  can  be  treated  with  the  same 
techniques. 
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The  mutual  coherence  function  (MCF)  of  the  acoustic  wavefunction  from  a  point  source  is  derived 
by  the  path-integral  technique  for  transmission  in  the  presence  of  a  sound  channel.  Separations  in 
time,  transverse  horizontal  position,  vertical  position,  and  acoustic  frequency  are  treated. 

Approximate  coherence  times,  lengths,  and  bandwidths  due  to  internal-wave  fluctuations  are 
derived.  The  MCF  of  frequency  is  explicitly  evaluated  for  fluctuations  due  to  internal  waves.  The 
shape  of  an  ensemble-averaged  pulse  is  derived  from  the  MCF  of  frequency. 

PACS  numbers:  43.30.Bp,  43.30.Ft,  43.20.Bi,  43.60.Cg 


INTRODUCTION 

The  mutual  coherence  function  (MCF)  contains  impor¬ 
tant  statistical  properties  of  the  acoustic  field  that  has  tra¬ 
versed  a  medium  filled  with  random  fluctuations.  For  exam¬ 
ple,  the  coherence  time  and  coherence  lengths  that 
determine  the  maximjim  effective  integration  times  and  ar¬ 
ray  lengths  that  can  be  utilized  in  sonar  systems  are  con¬ 
tained  in  the  MCF.  We  show  in  this  paper  how  to  derive  a 
number  of  results  about  the  MCF  which  were  indicated  in 
our  earlier  works  on  the  subject.1-2  In  particular,  we  show 
how  to  derive  the  MCF  from  the  path-integral  technique  for 
transmission  in  the  presence  of  a  deterministic  sound  chan¬ 
nel. 

The  MCF  of  frequency  controls  the  coherent  band¬ 
width  and  also  describes  the  behavior  of  the  ensemble-aver¬ 
aged  pulse  for  a  pulse-transmission  experiment.  In  this  pa¬ 
per,  we  derive  the  path-integral  expression  for  the  MCF  of 
frequency,  and  then  give  explicit  rules  for  calculating  it  in 
the  special  case  of  internal-wave  medium  fluctuations.  We 
are  able  to  explicitly  evaluate  the  case  of  internal  waves  be¬ 
cause  the  accepted  spectrum  of  internal  waves  implies  a 
nearly  quadratic  structure  function  and  path  integrals  for 
quadratic  actions  are  known. 

An  appropriate  Fourier  transform  of  the  MCF  is  the 
ensemble-averaged  pulse  (EAP)  at  the  receiver  for  a  narrow 
transmitted  pulse.  With  no  sound  channel  the  EAP  rises 
sharply  and  has  an  exponential  tail.  We  show  how  the  pres¬ 
ence  of  a  sound  channel  can  cause  a  precursor  region  to  the 
pulse,  and  can  cause  the  tail  at  late  time  to  decrease  in  inten¬ 
sity. 

I.  PATH-INTEGRAL  EXPRESSION  FOR  THE  ACOUSTIC 
WAVEFUNCTION 

We  begin  with  the  wave  equation  for  the  pressure  as  a 
function  of  space  and  time  in  the  presence  of  a  spatially  vary¬ 
ing  wave  speed.  We  follow  the  notation  of  the  review  article 
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by  Flatte3:  The  sound  speed  can  be  expressed  as 

C (x,t )  =  C0l  1  +  UJ[z)  +  n[x,t )) ,  (1) 

where  C„  is  a  reference  sound  speed,  UJz)  is  a  dimensionless 
function  of  the  depth  z  representing  the  deterministic  sound 
channel,  andjr(x,/ )  is  a  random,  zero-mean  function  of  posi¬ 
tion  representing  the  effect  of  medium  fluctuations  such  as 
internal  waves.  The  wave  equation  for  an  acoustic  wave  is 
unaffected  by  the  time  dependence  of  fs  because  fi  has  only 
components  with  very  low  frequency. 

The  acoustic  pressure  tp  obeys  the  wave  equation 

*»P-C2VV-  0.  (2) 

The  parabolic  approximation  consists  of  considering 
solutions  in  which  waves  arc  traveling  only  at  small  angles  to 
a  particular  direction;  in  the  ocean  this  direction  is  in  the 
horizontal,  labeled  by  x.  Thus  we  try 

P  =  exp[%*  -  ot )]  tfx,r  | ,  (3) 

where  q  and  a  are  the  wavenumber  and  frequency  of  an 
acoustic  wave  traveling  along  the  x  axis  at  speed  C0:  that  is, 
q  =  <r/Co.  The  “reduced”  wavefunction  is  slowly  varying 
in  space  and  time  compared  with  q  and  a,  and  satisfies  a 
parabolic  equation4-9 

+  =  [-<?„ -d*  +V(tfo+  A)J  *■  (4) 

Equation  (4)  is  a  Schrodinger  equation,  and  thus  its  so¬ 
lution  can  be  directly  expressed  in  terms  of  a  Feynman  path 
integral* 

iqSj(i)  -  iq  j£  h[x,j{x\ ,t  ]dx 

where  the  path  integration  (indicated  by  D )  is  over  all  paths 
4*1 =  ( yf*W*)]  connecting  the  source  to  the  receiver.  The 
phase  associated  with  the  path  in  the  absence  of  fluctuations 
is 

qS°  =  <l  £  l  +  ^  2)2  ~  0oW  ]  4x  ,  (61 

and  Nisa  normalization  factor  chosen  by  convention  so  that 
t/ma]  for /i  =  0. 
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„  path-integral  derivation  of  the  mcf-thi 
monochromatic  case 

The  MCF  for  the  smile-frequency  case  is  (  **(2)tHI)>. 
where  the  angle  brackets  indicate  averaging  over  the  ensem¬ 
ble  of  random  ft  functions.  The  MCF  measures  the  coher¬ 
ence  between  the  acoustic  fields  at  two  different  points.  la¬ 
beled  by  1  and  2.  These  points  may  be  separated  in  space, 
time,  or  both.  We  treat  spatial  separations  that  lie  in  the  ( yj) 
plane. 

The  use  of  (5)  and  (6)  results  in 

<  **(2M1» 

=  |2  J  Dz,  exp^fySoU)  -  iqSJfi) 

ft{z„ii)dx  +  ft{in,t1)dx^l .  (7) 

The  ensemble  average  applies  only  to  the/ts,  so  we  may  write 

<  ^*(2)  VKD> 

=  |J  J/fc,  Zfcj  exp{t*So(l)  -  iqSJ2)  -  J  Kl2) ,  (8) 
where 

yi2  =  q*([  jo  /a(x„r,)rfx  -  p^ytxj  ^  (9) 

and  we  have  used  the  fact  that 

<exp{/a)>  =  exp(  -  J<aJ)) ,  (10) 

if  a  is  a  zero-mean,  Gaussian  random  variable,  such  as  any 
combination  of  fis  is  assumed  to  be.  Even  if  a  is  not  a  GRV, 
( 10)  can  still  be  true  if  the  higher  moments  of  a  are  small. 

The  problem  of  finding  a  useful  result  for  the  MCF  now 
reduces  to  evaluating1  the  double  path  integral  in  (8).  We 
first  express  the  unperturbed  phases  in  terms  of  the  path 
variables  using  (6): 

So(l)  -  SJ2)  =  J[  [Jfc  J’.)1  -  «*,  ytf  +  J(<5,  *>)2 

-\(dxz7f-U^l)+U^2))dx,  (11) 

and  we  expand  the  deterministic  sound  channel  function  to 
second  order  in  the  displacement  of  the  paths  away  from  the 
equilibrium  ray.  That  is,  we  define  z,(x)  as  the  function  that 
satisfies  the  ray  equation  (in  the  parabolic  approximation): 


3„:,  +  UM,)  =  0  (12) 

with  boundary  conditions 

z,(0|  =  zt>  z,(*)=(z(l)-M2)]/2.  (13) 

Define  two  new  path  variables  by 

i>l(x)*z,(x)-zr(x),  VjIxJsszjM-z,!*),  (14) 

and  expand  (7o(z,)  around  the  point  z,(x): 

u&j  =  ujf,)  +  w0»x  +  i  uyt ,  (i5) 


where  it  is  understood  that  U J  and  I/q  are  evaluated  at  z,(x). 
This  expression  will  be  valid  as  long  as  the  effective  bundle  of 
acoustic  energy  stays  well  confined  around  the  unperturbed 
ray,  z,  (x).  If  more  than  one  solution  of  ( 1 2)  exists  (determinis¬ 
tic  multipath),  then  we  treat  one  unperturbed  ray  at  a  time. 
Addition  of  the  results  depends  on  the  coherence  between 


deterministic  rays,  which  can  often  be  well  estimated.  At 
caustics  our  method  breaks  down. 

In  terms  of  the  new  displacement  path  variables  till 
becomes 

5o(I)  -  5n(2|  =  Jf*  [J  13,  y,f  -  \[3,  y:f  +  Itf.u,)1 

-  +  3,zr\3,v,  -  3,v}) 

-  U'0{v,  —  w2)  —  |  U^v\  -  ii)]dx  .  Il6i 
Now  change  variables  once  again,  to 

a(x)  =  y,(x)  -  y2(x),  0{x)  =  (  y,[x)  +  y2(x))/ 2  , 
u(x)  =  u,(x)  -  u2(x)  =  z,(x)  -  z2(x) ,  |17| 

=  (t>,(x)  +  t»j(x)J/2  . 

Then  the  unperturbed  phase  term  becomes 

■Sod)  —  $o(2>  =  C  id,ad,  0  +  3,udxw  +  3X  z,3, u 

—  U gU  —  V o  uw]dx  .  (18) 

Now  the  crucial  observation  is  that  VK2  is  a  function  mainly 
of  the  difference  between  two  paths  z,  and  z2,  and  not  a 
function  of  the  average  (z,  +  z2)/2.  The  path  integral  (8)  is 
being  done  over  four  scalar  functions  a,  0,  u,  and  w.  The 
above  observations  correspond  to  noting  that  Vl2  is  not  a 
function  of  0  and  w. 

Consider  the  0  path  integral.  The  only  0  term  in  ( 1 8) 
gives 

I J  DP  exp^  (3,  ad,  0)dxj  .  ( 19) 

First,  integrate  by  parts  within  the  x  integral,  keeping  a  and 
0  fixed: 

Iff  =  j  DP  «p(*?  [  [  -  (<?„a)  P]dxj  .  (20) 

Now  the  P  integral  is  direct,  in  analogy  with  the  definition  of 
the  Dirac  delta  function,  and  yields  a  restriction  on  a: 

3„a  =  0.  (21) 

The  boundary  conditions  for  a  can  be  expressed  easily  for  a 
point  source: 

a(0)  =  0,  a[R )  =  y(l)  —  y{2) .  (221 

Thus  the  solution  for  a(x)  is 

a(x)  =  a{R  ix/R ) ,  (231 

which  is  the  separation  between  two  straight  lines  from  the 
source  to  the  two  receivers  at  positions  1  and  2.  Thus  the  a 
path  integral  is  effectively  done  despite  the  dependence  of 
V>2  on  a. 

Consider  the  w  path  integral.  First  we  can  integrate  the 
3,  z,dMu  term  by  parts  within  the  x  integral  yielding 
(  —  3„z,)u.  But  (12)  shows  that  this  cancels  the  U‘au  term. 
Integrating  the  3,  ud,  w  terms  by  parts  finally  yields 

So(l)  -So(2)  =  jM  -<?„«-  USu\u>dx,  (24) 

and  the  w  part  of  the  path  integral  yields,  again  in  analogy 
with  the  Dirac  6  function: 
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duu  +  u;u^  0.  1251 

The  boundary  conditions  for  u  can  be  expressed  in  the  same 
way  as  (22) 

«|0)  =  0.  utJl  )  =  r(l)-r)2).  |26) 

The  solution  of  (25)  and  126)  is  easily  shown  to  be  the  separa¬ 
tion  between  two  nearby  rays  in  the  sound  channel,  which 
start  at  the  source  and  end  at  the  two  receivers. 

Thus  the  path  integral  is  done,  and  the  only  require¬ 
ment  is  to  evaluate  P,2  at  the  separations  that  have  been 
determined  solely  from  the  integrations  involving  the  unper¬ 
turbed  phase  terms.  A  stationary-phase  approximation  has 
been  invoked  to  define  the  unperturbed  ray  in  the  absence  of 
fluctuations.  However,  the  fact  that  no  stationary-phase  ap¬ 
proximation  has  been  invoked  for  propagation  through  the 
fluctuations  should  be  emphasized.  But  K,2  for  z,(x)  and 
Zj(x)  being  nearby  rays  is  defined  as  the  well-known  phase 
structure  function: 

Vl2(a,u  =  ray  separation)==f>  (1,2).  (27) 

Thus  the  result  is 

<0*(2Wl)>=exp(-J/)(l,2)).  (28) 

At  this  point  it  is  worth  noting  that  a  number  of  end¬ 
point  terms  that  resulted  from  the  integrations  by  parts  have 
been  subsumed  in  the  normalization,  which  is  required  to 
give  unity  at  zero  separation  of  the  two  receivers. 

'  The  result  (28)  has  been  obtained  previously  in  many 
different  ways.*'11  The  path-integral  derivation  provides  a 
means  of  seeing  the  physics  of  the  approximations  in  a  new 
way.  The  only  approximations  have  been  the  parabolic  ap¬ 
proximation  and  that  Va  is  not  a  function  of  0  or  w.  The 
result  is  valid  in  the  presence  of  a  deterministic  sound  chan¬ 
nel,  regardless  of  whether  the  wave  fluctuations  are  unsatu¬ 
rated  or  not.  The  limitations  on  the  result  come  from  the 
approximation  being  violated  by  inhomogeneity  and  aniso¬ 
tropy.3 

III.  EVALUATION  OF  THE  MONOCHROMATIC  MCF  FOR 
INTERNAL  WAVES 


Under  the  Markov  approximation,  the  phase  structure 
function  can  be  expressed  as 


0(1,2)  =  Iq2  dx{n1)L,f{a,u,t), 


where  <  n 2)  and  Lf  have  been  defined  previously2  3  and 
f[a,u,t )  is  the  phase  correlation  function  (PCF)  defined  by 
Esswein  and  Flatte.12  The  PCF  has  been  evaluated  for  inter¬ 
na]  waves  by  Esswein  and  Flatte,12  using  a  combination  of 
analytical  and  numerical  techniques.  Since  a  and  u  are  both 
functions  of  range,  x,  (29)  must  in  general  be  evaluated  by  a 
numerical  integration  code. 

At  small  separations,  approximations  to  the  PCF  are 
possible.  For  example, 


/|0,0,/)  =  l|<u2!r2, 


(30) 


/(0,«,0)  =  J{k^(u2  \n\U(/u) ,  (31) 

where  |<t>2)  is  an  average  internal-wave  frequency  that  is 
dependent  on  the  local  depth  and  position  of  the  ray  zr{x), 
and  likewise  the  quantities  \k\.\  and  u0.  All  three  quantities 


are  evaluated  in  Esswein  and  Flatte.  Because  (31)  is  nearly 
quadratic  the  logarithm  may  rather  accurately  be  replaced 
by  a  constant.  The  best  constant  to  use  can  be  shown1 2  to  be 
In  <P.  where 

<P2  -  q2  jdx(p:)Lr  .  i32i 

Evaluations  of  <P  for  some  particular  examples  are  given  by 
Esswein  and  Flatte.13 

The  evaluation  of /fa, 0.0)  for  internal  waves  is  compli¬ 
cated  by  the  coupling  between  horizontal  internal-wave 
structure  and  both  vertical  structure  and  frequency.12  Nu¬ 
merical  evaluations  have  been  done12  and  seem  to  follow 


approximately  the  law 

/fa,0,0)s:  $  Bcf ,  (331 

where  the  power  p  is  empirically  about  1.5.  Thus  our  final 
approximate  results  for  the  monochromatic  MCF  are 

<  **[t  )<W)>  =  exp[  -  O.5(t/f0|2]  ,  (34| 

( tfi*(Az)  iMO))  =  exp(  -  0.5(dz/z0)2]  ,  135) 

<  i t>*\Ay)  $0)>  =  exp(  -  0.5(dy/yo)/’]  ,  (36) 

where  the  coherence  time  t „,  vertical  coherence  length  z0  and 
horizontal  coherence  length  y0  are  given  by 

‘o2=  r  J  dx(p-)Lr\a)'-\  ,  (37) 

Zo2=<f  lnV>f  dx(  /j.2)L,,\k\,  j  ,  (38) 

Jo 

dx(fi2)LrB,  (39) 


and  the  evaluations  of  all  the  quantities  in  (37H39)  can  be 
done  by  the  methods  of  Esswein  and  Flatte.'2 

Comparisons  of  (34)  and  (37)  with  experiment  are  done 
in  a  companion  paper.14  It  would  be  very  desirable  to  com¬ 
pare  (35)  and  (36)  with  experimental  data  in  an  appropriate 
parameter  range  in  which  internal  waves  should  dominate. 

IV.  PATH-INTEGRAL  DERIVATION  OF  THE  MCF  OF 
ACOUSTIC  FREQUENCY 

We  will  treat  in  detail  only  the  MCF  for  frequency  sepa¬ 
rations  with  no  simultaneous  space  or  time  separations. 
Then  (8)  is  replaced  by 

=  j,V  |2  J  Dt,  Dt,  expto.Sol  1)  -  *A(2)  -  J  Ki; ] . 

401 

where  q,  =  <x,/C0  and 

V'2  =  ((?I  i  VdvWx  -  q2  /i(z:,0)<fr)2J  .  (41i 

The  unperturbed  phase  is  expressed  in  the  manner  of 
(1 1)  as 

*AIU  -  9.'5o(2)  =  [1  qx(dxyx)2  -  l  ?-(<?, y-r 

-  19i!«9,  z,)2  -  1  ?:(<?,  Z;): 

-^tfoM  +  fjUolZjJldx.  1421 
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Applying  the  same  changes  of  variables  indicated  in  ( 14) 
and  (17),  we  And 

i)  -  «Ai2)  -  ?i*n>  -  sm  +  w.0)1 

+  i ft, a)2  +  (d.  ur)2  +  l(dMuf 

+  2dMwdt2r  +  (dMz.)2-2U^r) 

-  2U’0w  -  VW  -\UZu2]dx,  (43) 

where  q  =  ( q,  +  4j)/2,  Aq  =  qt  —  4*  and  SJl)  -  S„( 2)  is 
expressed  in  these  variables  in  (18).  Again  we  observe  that 
Ki2  is  not  a  function  of  0  and  w,  and  this  allows  the  path 
integrals  over  0  and  u>  to  be  done.  However,  in  this  unequal 
frequency  case  the  path  integrals  are  not  analogous  to  6  func¬ 
tions;  instead  they  are  Gaussian  integrals  because  of  the  qua¬ 
dratic  forms  in  (43). 

Consider  the  0  integral.  Collecting  0  terms  in  (43)  we 

have 

Ifi=  J  D0  e*p[y  44  J[  (&  P  )J  +  ^  d*ad*  &  )  <**]  • 

(44) 

This  path  integral  is  done  by  completing  the  square,  and, 
aside  from  factors  that  will  be  subsumed  in  the  normaliza¬ 
tion,  yields 

Now  consider  the  w  path  integral.  This  can  be  expressed 
as 

/„  *  J  Dw  exp|y  [wLw  +  ^  j  . 

(46) 


where  L  is  an  operator  given  by 

L=-d„-U;.  (47) 

Note  that  integration  by  parts  has  been  used  in  several  places 
to  obtain  (46).  Evaluation  of  /„  is  done  by  completing  the 
square  and  yields,  apart  from  normalization  terms, 

A--K*r~4  1481 

Combining  (43),  (43),  and  (48),  we  find  (40)  becomes 
<  **(  qiM  4i)>  =  AT  J  flo  Z>w  exp(  -  y 

xj[  [(dxa)2  +  uLu)dx  -  yFl2). 

(49) 

First,  let  us  note  that  if  Aq  =  0,  then  this  path  integral  is 
essentially  a  stationary-phase  integral  around  the  extremum 
of  the  range  integral  in  the  exponential.  The  solutions  for  a 
and  u  are  then  (23)  and  (25)  and  (26),  respectively.  This  MCF 
expression  (49)  is  a  path  integral  over  two  scalar  (or  one  vec¬ 
tor)  path,  and  it  cannot  be  simplified  without  some  assump¬ 
tions  about  the  nature  of  K,2.  Let  us  express  Fl2  in  terms  of? 
and  Aq: 


k,j = **({  I -  j 

+  4Jf[((  J MfciKfr)  )  -  ((  J ^(*jkfr)  }] 

+  y  lAqf^  j  fi{tt)dx  +  J A*(*2)d*)  }  •  (501 


The  lack  of  dependence  of  K,2  on  the  centroid  of  the  paths 
( 0,w)  allows  us  to  neglect  the  qAq  term  and  to  estimate  the 
last  term  as  (Aq/  q)2<t> 1 ,  where 


*i=  **((  J )  ' 


(51) 


Since  </> 1  is  a  function  only  of  the  unperturbed  ray  i„ 
and  not  the  paths,  the  path-dependent  part  of  K,2  is  only  the 
first  term  of  (30).  Using  the  Markov  approximation  at  this 
point  we  have  (32)  and 

Vl2  =  (Aq/qfd>2  +  2q1  f  {n2)Lrf(a,u,0)dx.  (52) 


The  dependence  on  paths  is  contained  in/(a,u,0),  and  the 
expression  for  the  MCF  is  now 


< 


^*(4jM4 i)>  =  exp  -  y  * 2]  Q  M4) »  l53> 


Q(Aq) 


=  ivj  AxDuexp^- yy~J[  [(dxa)2  +  ul.u]dx 
-rf  {n2)Lrf(a,u,O)0xj. 


(54) 


The  evaluation  of  the  path  integral,  Q  (Aq),  depends  on 
our  understanding  of  f(a,u,0).  We  make  use  of  the  approxi¬ 
mations  given  in  (31)  and  (33),  and  our  knowledge  of  the 
magnitudes  of  B  and  |k  \  |,  to  justify  the  neglect  of  the  de¬ 
pendence  of/on  the  horizontal  variable  a  (because  B  is  very 
small).  The  a  integral  then  yields 


Q(Aq)  =  N 


uLudx 

Y^ln^jf  (n2)Lr\klv\u2dxj. 


(55) 


This  is  a  one-dimensional  path  integral  with  quadratic 
action,  whose  explicit  evaluation  will  be  given  in  the  next 
section.  Note  that  the  normalization  N  is  set  so  that 
Q  (Aq)  =  l  if fi  —  0  everywhere. 


V.  EVALUATION  OF  THE  MCF  OF  FREQUENCY 

The  evaluation  of  Q  (Aq)  from  the  path  integral  (55)  can 
be  done  in  at  least  two  useful  ways.  The  first  involves  fixing 
Aq  and  solving  an  ordinary  differential  equation  for  Q.  The 
second  involves  an  eigenvalue  method  that  finds  the  singu¬ 
larities  in  Q  (Aq). 

It  will  be  useful  to  express  Q  in  the  following  way: 

Du  exp^  - 

Af(x)=  -du-UZ-iAq\n<P{n2)Lt\k2v].  (57) 


-Lq2^  uMu  dx^  ,  (561 
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where  we  have  used  the  feet  that  s  q1  ■  The  solution  of 
thi«  quadratic  path  integral  is  given  in  terms  of  the  solution 
of  a  differential  equation*  for  a  function  S  \x^q).  The  equa¬ 
tion  is 

MS  =  0,  (58) 

with  initial  condition  S  ((Xdq)  =  0.  Then 

Q[Aq)  =  [S\R,0\/S(R4q))ul  ■  (59) 

Note  that  5  is  complex. 

The  form  (59)  is  useful  for  finding  the  behavior  of  the 
MCF  at  small  frequency  separations.  (Remember  that  the 
frequency  separation  Aa  =  AqC0,  and  we  define  qC0.) 
Consider  M  as  made  up  of  two  parts 

M  =  L  —  AaM, ,  (60) 

Af,  =  iC0~  1  In  ♦  ( fi2)LF\ky  \  .  (61) 

Considered  in  position  space  these  operators  are  infi¬ 
nite-dimensional  matrices,  and  the  inverse  of  L  can  be  found 
by  the  equation 

(  -  -  US)L  ~  W')  =  Six'  -  x) ,  (62) 

hence, 

L~\x*‘  )=«(x^'),  (63) 

where  g{xjc’)  is  the  Green’s  function  defined  in  Ref.  3;  it 
depends  on  the  sound  channel  through  U 

At  small  Aa,  we  may  solve  (58)  by  a  perturbation  expan¬ 
sion.  This  yields 

ta(e(d<r)]=^Tr<Z.-,A/I) 

+  ^£}L  h[L-'M,L-'Sfx),  (64) 

4 

where  Tr  indicates  the  trace.  The  traces  become  integrals  in 
the  infinite-dimensional  case,  so  that  we  may  finally  express 
Q  as 

ln[Q(do)]  ssidor,  -  J(d<rr0)2  ,  (65) 

IH!  fc*\4x,  (66) 

-fe) 

X  dx'lSyLrlkUlglxS)]1.  (67) 

where  it  is  understood  that  the  first  ( fi2)L,  ( k  {, )  is  evaluat¬ 
ed  at  x  and  the  second  at  x‘. 

Although  (65)  is  only  exact  as  Aa— *0,  it  is  worthwhile  to 
model  the  MCF  as 

(  =  exp |  -  y  [~J  <P 2  j 

X  exp(/d<7T,  —  }(dcrr0)2j  ,  (68) 

where  the  second  exponential  comes  from  the  Q  function.  In 
this  form  it  is  clear  that  r,  is  a  shift  of  the  time  origin  for  the 
frequency  determination,  and  r0is  to  be  compared  with  <t>  /e r 
to  find  which  term  is  more  important  in  determining  the 
width  of  the  MCF.  That  is,  r0" 1  and  o/<P  combine  to  deter¬ 
mine  the  coherent  bandwidth  of  the  transmission. 


TABLE  I.  Parameter!  of  QtAa I  for  (wo  example*  with  constant  F Ixi  and 
Un  full*  set  at  10  1  msAiti1  Unit*  lor  the  rs  art  ms  Untu  for  the  eigen¬ 
value*  arc  im»i " 1  The  eigenvalue*  were  calculated  by  a  general  computer 
code  (CO;  they  should  agree  with  <721  for  these  ample  examples,  which  they 
do  reasonably  well. 


Uotkm  ' 

0 

0.034 

r, 

1.02 

-2.61 

U 

0.65 

2.76 

cc 

(72) 

cc 
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4. 

0.80 

0.81 

-2.60 

-2.59 

3.22 

3.24 

—  018 

-0.18 

7.24 

7.29 

3.84 

3.85 

12.87 

12.89 

9.47 

9.49 

^5 

20.1 

20.1 

16.71 

16.47 

Ajo 

80.4 

80.6 

77.0 

772 

^15 

181 

181 

178 

178 

^20 

322 

322 

318 

319 

The  behavior  of  the  MCF  near  zero  Aa  is  only  part  of 
the  information  available.  To  learn  more,  we  return  to  (58) 
and  solve  it  by  an  eigenvalue  method: 

LS.(x)  =  /l.F(x)S.W.  (69) 

/»  =  Co'1  In  *  </r2>Z,,  |*i.  |  ,  (70) 

with  boundary  conditions  S„(0)  =  S„(R  )  =  0.  Then 

QlAa)  =  \\[\~iAa/Xn\-m  ■  (71) 

A 

There  are  many  numerical  methods  for  solving  (69).  We 
find  the  most  effective  to  be  expanding  Sm  (x)  in  Fourier 
modes  along  with  the  eigenvalues  4. .  In  applying  this  meth¬ 
od  one  must  be  sure  that  enough  Fourier  modes  to  accurate¬ 
ly  represent  the  S„  up  to  the  desired  maximum  n  have  been 
taken. 

Both  F(x)  and  U  £  affect  the  values  of  A. .  As  an  exam¬ 
ple,  take  U  o  and  F[x )  =  constant.  Then 

X,  =F-'((jr/*)V  -US),  n  =  1,2, .  (72) 

Note  that  all  of  the  A.  are  positive  if  U£  =  0  (no  sound  chan¬ 
nel).  In  general,  a  nonzero  U  J  can  pull  the  lowest  lying  X, 
below  zero.  The  number  of  X .  below  zero  is  equal  to  the 


FIG.  1.  The  MCF  of  frequency  for 
no  sound  channel  and  constant 
F(x)  =  10" 1  tmsl/km1.  The  range  is 
taken  to  be  35  km.  The  solid  cum  is 
calculated  by  a  general  computer 
code.  The  dotted  curve  is  the  approx¬ 
imation  |6SI  with  |66|  and  (67),  which 
give  r,  =*  1.02  ms  and  r0  =  0.65  ms 
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FIG.  2.  The  same  as  Fig.  I,  except 
the  sound  channel  is  represented  by 
constant  UZ  »  0.034  (kml'1. 

r,  *  —  2.61  ms.  and  r0  »  2.76  ms. 


number  of  caustics  the  unperturbed  ray  has  passed  from 
source  to  receiver.15  As  a  particular  /l„  passes  through  zero, 
the  receiver  is  passing  through  a  caustic. 

Consider  two  examples  in  which  4>/a  =  2.8  ms, 
F—  10-i  ms/km2,  and  either  =  0,  or  0.034  km-2. 
These  values  are  comparable  in  magnitude  to  those  in  the 
AFAR  experiment14  which,  however,  has  an  Fand  U  J  that 
vary  with  range.  The  range  is  taken  to  be  35  km.  Table  1  gives 
the  values  of  the  eigenvalues  for  the  examples,  as  calculated 
from  our  numerical  solutions  to  (69).  They  agree  with  the 
analytic  expectation  (72).  (The  same  computer  code  has  been 
used  in  calculating  the  eigenvalues  for  the  AFAR  experi- 
ment. M)  Figure  1  shows  the  Q  [Aa)  for  this  case.  The  values  of 
r,  and  t0  calculated  from  (66)  and  (67)  are  given  in  Table  I. 
These  values  (it  the  phase  slope  and  the  l/e1/2  magnitude 
half-width  rather  well. 

When  U  £  is  set  to  0.034  km  ~ 2,  the  resulting  Q  ( Aa )  has 
quite  a  different  phase  behavior;  the  phase  slope  has  changed 
sign  (Fig.  2).  The  significance  of  this  slope  is  discussed  in  the 
next  section.  The  values  of  r,  and  r0  calculated  from  (68)  and 
(67)  are  given  in  Table  I,  for  this  case,  and  again,  they  fit  the 
phase  slope  and  magnitude  half-width  rather  well. 

VI.  PULSE  TRANSMISSION 

Consider  the  Fourier  transform  of  the  MCF  as  a  func¬ 
tion  of  4<r,  label  the  conjugate  variable  as  r,  which  we  call 
pulse  time.  If  an  experiment  sends  many  pulses  from  a  fixed 
source  to  a  fixed  receiver,  then  the  shape  of  the  ensemble- 
averaged  pulse  (EAP)  is  the  Fourier  transform  of  the  MCF. 

What  are  the  meanings  of  the  various  parameters  we 
have  used  to  characterize  the  MCF,  such  as  r ,,  and  the  set 
of  4.  ?  First,  r ,  is  a  shift  in  the  mean  arrival  time  of  the  EAP. 
That  is,  the  medium  fluctuations  cause  a  mean  shift  of  the 
pulse  from  when  it  would  arrive  in  the  absence  of  fluctu¬ 
ations,  and  that  shift  is  r,.  Unfortunately,  very  few  experi¬ 
ments  (none  in  the  ocean)  can  measure  the  pulse  time  in  the 
absence  of  fluctuations,  and  the  absolute  sound  speed  is  not 
known  well  enough  to  calculate  it,  so  a  prediction  of  r,  can¬ 
not  t>.  directly  tested.  Second,  r0  is  a  contribution  to  the 
avuag.  widt1-  of  the  pulse,  as  is  <P  /a.  Both  act  as  standard 
deviations  in  a  Gaussian  EAP. 


The  tails  of  the  EAP  at  both  early  and  late  times  ?re 
controlled  by  the  smallest  X . .  This  can  be  seen  by  noting  that 
(7 1 )  is  a  product  of  terms,  so  its  Fourier  transform  is  a  convo¬ 
lution  of  the  Fourier  transforms  of  all  the  terms  separately. 
In  a  convolution,  the  tat!  is  controlled  by  the  term  with  the 
widest  Fourier  transform.  Therefore  we  need  the  Fourier 
transform  of 

&,(4<r|  =  (l  -LAa/XJ-1'2.  (73) 

where^m  is  the  eigenvalue  with  the  smallest  magnitude  with 
the  appropriate  sign.  The  Fourier  transform  is 


PJr)  =  {Xm/irT)‘lle~'imr,  Xmr>  0, 
=  0,  Xm  r< 0, 


(74) 


where  we  have  normalized  so  that  /  ( I{r))dr  =  1.  To  find 
the  shape  of  the  tail  we  convolve  this  function  with  the  rest  of 
the  MCF:  that  is,  with  a  Gaussian  of  standard  deviation  <P  / 
a ,  and  all  the  other  eigenvalue  terms,  resulting  in  a  tail  given 
by 


— [l(¥)][n(-^r] 


(3  \  1/2 

-=-]  e-A-T,  Xm  r>0 .  (75) 


The  smallest  positive  eigenvalue  dominates  the  very  late  r 
behavior.  The  smallest  negative  eigenvalue  gives  a  precursor 
and  dominates  the  very  early  r  behavior.  If  there  are  no  nega¬ 
tive  eigenvalues  (a  caustic  has  not  been  passed)  there  is  no 
precursor. 

Consider  the  examples  given  in  Figs.  1  and  2.  The  corre¬ 
sponding  EAPs  are  shown  in  Figs.  3  and  4  along  with  the 
predicted  tail  and  precursor  from  (75).  Note  that  the  approx¬ 
imation  for  Q  given  by  (68)  would  give  the  position  and  gross 
width  of  the  pulse,  but  would  not  give  the  precursor  in  Fig.  4. 
Note  also  that  if  the  controlling  eigenvalue  has  a  large  mag¬ 
nitude  (such  as  X3  =  +  3.85  ms  ~ 1  for  the  positive  r  tail)  then 
the  Gaussian  remains  in  control  to  much  larger  r.  The  time 
at  which  (75)  would  become  relevant  is  approximately 
rs;Xm  (<P  /a)2,  or  30  ms  in  the  case  of  Fig.  4. 

The  EAP  calculation  is  valid  whether  the  transmission 
is  unsaturated  or  not.  In  the  unsaturated  region  a  narrow 
transmitted  pulse  will  be  received  as  narrow,  and  the  precur¬ 
sor  or  tail  of  the  EAP  is  formed  by  there  being  a  probability 


1  JHT-r— -Hr— — TTTTTTTTT 
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FIG.  3.  Ensemble-avenged  pulse  calculated  for  the  example  given  in  Fig.  I . 
The  solid  curve  is  calculated  by  a  general  code.  The  dotted  curve  is  from 
|68|.  The  dashed  curve  is  the  prediction  of  the  tail  from  (73)  where  the  con¬ 
trolling  eigenvalue  is  whose  value  is 0.81  ms". 
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FIG.  4.  Ensemble-averaged  pulse  calculated  for  the  example  given  m  Fig.  2. 
The  curves  have  the  same  significance  as  in  Fig.  3.  The  controlling  eigenval¬ 
ue  is /t,  (for  negative  r'  hose  value  is  —018  ms”'. 

for  the  pulse  to  be  early  or  late,  as  well  as  a  probability  for 
pulses  arriving  at  different  times  to  have  different  intensi¬ 
ties.5  In  the  saturated  region  each  pulse  may  have  its  individ¬ 
ual  precursor  or  tail. 

VII.  SUMMARY  AND  CONCLUSION 

The  derivation  of  the  mutual  coherence  function  (MCF) 
of  time,  transverse  space,  and  frequency  by  the  path-integral 
technique  has  been  given.  Allowance  for  a  deterministic 
sound  channel  and  the  presence  of  reasonable  inhomogene¬ 
ity  and  anisotropy  in  the  fluctuation  field  has  been  included. 
The  MCF  has  been  evaluated  for  fluctuations  dominated  by 
internal  waves,  which  have  a  vertical  structure  function  that 
is  nearly  quadratic.  Reasonably  simple  expressions  in  terms 
of  environmental  measurements  for  acoustic  coherence 
times  and  lengths,  and  coherent  bandwidths  have  been  giv¬ 
en. 

The  Fourier  transform  of  the  MCF  of  frequency  is  the 
ensemble-averaged  pulse  (EAP).  The  possibility  of  a  precur¬ 
sor  or  a  tail  on  the  EAP  is  shown  to  depend  on  both  the 
deterministic  sound  channel  and  the  fluctuations,  and  meth¬ 
ods  have  been  given  to  calculate  these  effects. 


All  of  the  parameters  that  can  be  used  to  approximate 
the  various  MCFs  and  the  EAP  can  be  evaluated  on  a  VAX 
computer  in  a  few  seconds  for  a  typical  ocean-acoustic  ex¬ 
periment. 
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Abstract 

A  series  expression  is  developed  for  the  fourth  moment  of  a  beamed  field 
incident  on  a  random  phase  screen  or  an  extended  medium.  The  series  has  a 
symmetry  that  allows  its  first  few  terms  to  generate  useful  approximations  at 
both  low  and  high  spatial  frequency.  The  parabolic  wave  equation,  the  Markov 
approximation,  and  Gaussian  refractive  index  fluctuations  are  assumed.  The 
result  for  the  phase  screen  is  obtained  by  Green’s-function  techniques.  The 
extended-medium  result  is  derived  in  an  analogous  manner  using  path  integral 
methods.  The  same  results  are  also  derived  by  moment-equation  methods.  The 
behavior  of  the  leading  terms  is  compared  to  previous  results  for  plane-wave  and 
point-source  geometries. 
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1.  INTRODUCTION 


Many  years  of  research  have  been  devoled  to  the  study  of  wave  propagation  in  random 
media  (WPRM).  The  first  comprehensive  review  of  the  field  was  by  Tatarskii  [1971],  followed  by 
Prokhorov  et  al  [1975],  Ishimaru  [1978],  and  Fante  [1975,1980].  The  propagation  of  radio 
waves  through  the  ionosphere  is  reviewed  by  Yeh  and  Liu  [1983].  The  phenomenon  of  interstel- 
lar  scintillation  is  reviewed  by  Rickett  [1977]  and  Rickett  et  al  [1984],  Sound  propagation 
through  the  ocean  and  path  integral  techniques  are  discussed  by  Flatte  et  al  [1979]  and  Flatte 
[1983]. 

We  consider  waves  propagating  from  an  arbitrary  source  distribution  in  a  random 
medium.  We  assume  the  statistics  of  the  medium  are  locally  homogeneous,  and  we  make  the 
Markov  approximation;  i.e.  the  field  fluctuations  induced  within  a  correlation  length  along  the 
propagation  direction  are  weak.  For  a  more  complete  discussion  see  Codona  et  al,  [1985]. 
The  wave  propagation  is  characterized  by  narrow  angular  scattering  due  to  the  small  random 
fluctuations  in  refractive  index.  It  is  then  convenient  to  write  the  complex  monochromatic 
scalar  field  as  E{%z)eika  where  *  is  the  propagation  direction,  iis  the  transverse  coordinate 
and  k  is  the  wavenumber  of  the  wave  with  no  refractive  index  fluctuations. 

The  random  nature  of  the  fields  is  conveniently  described  by  statistical  moments 
evaluated  in  the  transverse  plane  located  at  distance  /?.  Ensemble  averages  of  random  vari¬ 
ables  are  denoted  by  <>.  The  first  moment 

T  ,<**)  =  <£•(*.*)>  (1) 

or  average  of  the  field  and  the  second  moment 

)  =  <E(tx,R)l?  (**.*)>  (3) 

or  mutual  coherence  function  are  well  understood  [Tatarskii,  1971].  However,  there  are  few 
analytic  results  for  the  fourth  moment 

r4(ll,te,la.*4./?)  =  <E{tx.R)tf  fa,R)E(t3,R)£  (*«.*)>  (3) 

Previous  theoretical  work  concentrated  on  plane-wave  and  point-source  geometry.  We 


present  three  main  results  for  arbitrary  source  distribution. 

A  series  expression  for  the  fourth  moment  is  derived  as  an  expansion  of  the  Green's  func¬ 
tion  for  the  fourth  moment,  thus  avoiding  the  difficulties  associated  with  the  source  distribu¬ 
tion.  For  the  thin-screen  problem,  the  expansion  quantity  is  a  combination  of  phase  structure 
functions.  For  the  extended  random  media,  the  expansion  quantity  is  an  analogous  combina¬ 
tion  of  phase  structure  function  densities.  The  Green's  function  is  expressed  as  a  multiple 
path  integral.  The  resulting  series  of  path  integrals  is  evaluated  with  a  useful  identity. 

Our  second  result  is  the  generation  of  two  series  for  the  intensity  correlation  or  intensity 
spectrum.  The  fourth  moment  has  the  obvious  symmetries  that  it  is 

unchanged  by  interchanging  and  ia  or  by  interchanging  3^  and  i|.  Each  term  of  the  series 
expansion  does  not  share  the  symmetry  of  the  entire  expression.  Thus  two  separate  series 
are  obtained  by  invoking  symmetry.  In  principle,  either  series  could  be  summed  to  give  1%.  We 
demonstrate,  however,  that  it  is  better  to  consider  both  series  in  order  to  describe  the  fourth 
moment  with  the  fewest  number  of  terms.  This  assertion  is  demonstrated  for  the  second 
moment  of  intensity  or  intensity  correlation,  C(it,$2,R).  \'hich  is  a  special  case  of  the  fourth 
moment,  i.e. 

c(*, .**/?) = </(tl,f?)/(*z,/?)> = r4(tI,i1,i2.i2) = r4(ii.*.t M  (4) 

Note  that  the  symmetry  of  the  fourth  moment  has  been  explicitly  indicated.  A  clear  presenta¬ 
tion  of  the  behavior  of  the  intensity  correlation  series  obtained  from  the  fourth  moment 
expansion  requires  the  introduction  of  a  spatial  spectrum  of  intensity  fluctuations  for  a  spa¬ 
tially  nonstationary  random  process.  We  adopt  the  definition 

*<AA*)  =  dp  (5) 

where 

A=  ~(*i  +  <e)  P  =  * ,-tg  (6) 

(Note  the  free  format  of  the  argument  list  of  functions).  The  spectrum  has  the  property 

m 

f  *(JLA'R)  d^=C{ft.O.R)  =  <l(ft,R)*>  (7) 

- m 

It  should  be  noted  that  the  spatial  spectrum  may  depend  on  the  centroid  ft. 
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Since  there  ere  two  series  for  the  intensity  correlation  there  are  also  two  series  for  the 
intensity  spectrum.  The  leading  terms  of  one  series  for  describe  the  small  )) 

behavior  while  the  other  series  is  valid  at  high  Q.  The  rate  of  convergence  of  each  series  pro* 
vides  a  criterion  for  merging  the  two  results  to  produce  a  complete  expression  for  the  inten* 
sity  spectrum.  In  general,  an  analogous  treatment  of  the  intensity  correlation  series  is  not 
possible  since  the  leading  terms  of  both  series  do  not  converge  to  the  variance  as  the  spatial 
separation  approaches  zero. 

Our  third  result  is  the  demonstration  of  the  equivalence  of  path  integral  and  moment- 
equation  methods.  Early  theoretical  work  on  WPRM  concentrated  on  geometrical  optics  and 
the  method  of  small  perturbations  [Barabanenkov,  1971;  Tatarskii,  1971].  These  two 
approaches  were  limited  to  weak  scattering  conditions.  This  restriction  was  removed  with  the 
introduction  of  differential  equations  for  the  moments  of  the  field  [Prokhorov.  1975].  Func¬ 
tional  techniques  of  high  energy  physics  (path  integrals  and  operator  methods)  provided 
another  point  of  view  to  WPRM  [Klyatskin,  1973;  Dashen,  1979].  The  moment  equation  method 
and  functional  techniques  are  equivalent  [Codona  et  al,  1965]  and  must  generate  identical 
results  when  expansions  are  performed  in  the  same  quantity.  This  equivalence  is  demon¬ 
strated  by  deriving  the  same  fourth  moment  series  expression  using  moment-equation 
methods. 

The  thin-screen  case  is  considered  in  section  2.  The  second  and  fourth  moment  are 
analyzed  with  Green's  function  techniques  and  the  behavior  of  the  intensity  correlation  is 
investigated.  The  same  analysis  for  the  extended  medium  case  is  presented  in  section  3.  Here 
we  use  the  path  integral  representation  of  the  Green's  function.  Identical  results  for  the 
fourth  moment  are  derived  with  moment  equation  methods  in  section  4.  The  main  results  of 
the  paper  are  summarized  in  section  5. 
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Z.  GREEN  S  FUNCTION  APPROACH  TO  THE  THIN  SCREW  PROBLEM 


2.1  Introduction 

One  of  the  first  WPRM  problems  considered  was  the  propagation  of  plane  waves  through  a 
random  phase  screen  [Mercier,  1962;  Salpeter,  1967;  Bramley,  1967;  Gochelashvily  and 
Shishov,  1971,  1972,  1975;  Rumsey,  1975;  Rino.  1979a,  b;  Uscinski  and  Macaskill,  1963a.  b]. 
The  propagation  of  radio  waves  through  the  ionosphere  and  the  solar  wind  are  two  applica¬ 
tions  of  this  model.  The  theory  of  scintillation  from  a  point  source  viewed  through  a  random 
phase  screen  has  been  investigated  by  Lee  [1977].  The  case  of  a  Gaussian  beam  focussed  on 
the  observation  plane  has  been  considered  by  Gochelashvily  [1974].  Previous  work  concen¬ 
trated  on  plane-wave  and  point-source  geometries.  We  analyze  the  more  general  problem  of 
an  arbitrary  beam  incident  on  a  phase  screen  using  Green's  function  methods.  The  following 
analysis  is  presented  in  a  fashion  that  permits  a  clear  extension  to  the  more  complex  problem 
of  wave  propagation  in  extended  random  media.  We  review  Green’s  function  methods  with  a 
discussion  of  the  second  moment.  A  series  expression  for  the  fourth  moment  is  presented  as 
an  expansion  of  the  Green's  function  for  the  fourth  moment.  The  behavior  of  the  resulting 
series  for  the  intensity  correlation  is  then  discussed. 

Consider  the  scalar  wave  field,  E( f,*).  incident  on  a  thin  random  phase  screen  situated 
at  the  plane  z=0.  The  field,  E{  f.O  +  ),  emerging  from  this  screen  is  given  by 

^f.O*  )  =  £•(*,  0)e<e(t>  (8) 

where  E(t, 0)  is  the  field  just  before  the  interface  of  the  screen  and  the  phase  fluctuations  are 

o» 

&(r)  =  kfn(r.z)dz  (9) 

0 

where  n(?,z)  is  the  random  fluctuations  in  refractive  index.  Assume  8(i*)  is  a  zero  mean 
Gaussian  random  variable  with  homogeneous  statistics  and  correlation  function 


c*(t)  =  <e<fr)e{*  ♦  i)> = d* 
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where  *#(4)  is  the  spectrum  of  phase  fluctuations.  The  structure  function  of  phase  fluctua* 
tions,  #•(£)•  is  defined  by 

m 

D*(fi)  =  <[0(t)-6(t +  ?)]*>  *2[C.(0)-C.(^))  =  2/[1  -cosft.  £)]♦,(*)  rf*  (11) 

For  narrow  angular  scattering  the  scalar  field  satisfies  the  parabolic  wave  equation 

2i*||-+Vz£’  =  0  (12) 

The  solution  of  the  field  at  a  distance  R  from  the  screen  is  (Mercier,  1962) 

«» 

E(%R)  =  f  E(t,0)G(tt.R)  dt  (13) 

where  the  Green’s  function  is 

G(t*.R)  =  e*WGf(i;*.R)  (14) 

with  the  free  space  Green’s  function 

G'(*r*)=—  (15) 


2.2  Green’ s-F\inction  Approach  to  the  Second  Moment 

Since  the  random  fields  have  a  Green's  function  solution,  the  moments  of  the  field  also 
have  a  Green’s  function  representation.  Consider  the  second  moment  of  the  field  or  mutual 
coherence  function.  JVjfj.i*,./?),  given  by 

•»  m 

rs{tl.t2.R)  =  <E(ti,R)£  {te,R)>~  ffrd*i.r2,0)<G(ll,rx,R)(?(te;*i.R)>dtl  dt2  (18) 

The  Green's  function  for  the  second  moment  is  identified  as  the  expression  inside  the  <>  ,  i.e. 

G2(ii;ft.f?)  =  G{(it;f'1./?)<exp{i(d(i'l)-d(i'2)))>  (17) 

where 

Ci(*i-.*t.R)  =  & (!,;*,.*)</  {*t.tt.R)=  jg^exp[|^[(i,-i'Of -(*-*«)*])  (18) 
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is  the  free  space  Green's  function  for  the  second  moment.  Here,  t%  and  It  denote  the  set  of  f 
and  t  coordinates  respectively.  The  expectation  over  the  random  phase  is  performed  with  the 
identity 

-f  <e*> 

<e<#>  =  c  •  (19) 

which  is  valid  for  zero  mean  Gaussian  random  variables,  The  Green's  function  for  the  second 
moment  is  then 

G^rt.*^  Gf(*i;*'t,7?)exp[-^l?e(*'i-*'2)]  (20) 

2.3  Green* s-f unction  Approach  to  the  Fourth  Moment 

Previous  work  on  the  fourth  moment  has  concentrated  on  plane-wave  conditions 
[Zavorotnyi,  1979a].  We  address  the  full  fourth  moment  with  an  arbitrary  source  distribution. 
Following  the  previous  analysis  of  the  second  moment,  the  Green's  function  for  the  fourth 
moment  is  given  by 

G4(ii;fi./f)=G((il;ir4.f?)  (21) 

e*p[  -  “[£•(*,  -  fa)  +  0a(f  s  -  *«)  +  0#(f ,  -  f  ♦)  +  D9(tt  -  f  3)  -  D9(t2  -  f  4)  -  De(t ,  -  *,)]] 

where  G{ is  the  free  space  Green's  function  for  the  fourth  moment.  It  is  convenient 
to  apply  the  unitary  coordinate  transformation  [Rumsey,  1975] 

2d  =  i,  +  *g  +  *a+i«  2ii  =  d  +  £  +  t  +  $ 

2^1,  +  ^-13-i,  2*8  =  a  +  j5-f-i 

2?  =  i,-*g-ia+»4  2f8  =  d-*-?  +  « 

23  =  »,-tg  +  i3-*4  2*4  =  a-$+?-3  (22) 


The  set  (d'.^'.f ,3')  ,  will  denote  the  same  transformation  on  the  coordinate  set  (ft.f2.f3.f4). 
The  Green's  function  for  the  fourth  moment  then  becomes 
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g&;Tc.*)-g(  (*;*«.*)  (23) 

exp  - ^[Pi(r  ♦  3#)  ♦  dm  -i)  ♦  £•(*•  ♦  *)  ♦  Dtf  - «•)  - Dtf  ♦  r>  - />•(*•  -r » 

where 

(K(ti:t,.R)  =  -*'i),-(«e  -  **)*  ♦  «,-*,)*-(*. -**«)*) 

=  p^;<!xpU-I(a-a').(2-2')t(^-^).«-r))  (Z4) 

is  the  free  space  Green's  function  for  the  fourth  moment.  This  expression  is  intractable,  both 
analytically  and  numerically.  Marians  [1975],  numerically  calculated  the  intensity  spectrum 
of  plane-waves  incident  on  a  two-dimensional  phase  screen. 

The  plane-wave  case  was  considered  by  Zavorotnyi  et  al  [1977].  They  noted  that  a  combi¬ 
nation  of  structure  functions  was  small  in  the  important  regions  of  integration.  A  series 
expression  was  then  obtained  as  a  Taylor  series  expansion  For  the  general  case,  we  identify 
that  same  expansion  quantity  as 

e=  -j[D*(fi‘  +  l)  +  D9&-l)-D9@  +  r)-Dt&-r)]  (85) 

=  2  f  *e(*)[cos{  J ’  •  t)  -  cos(?  •  2)  Je4*  **  d  Z 


We  will  investigate  the  behavior  of  Q  for  a  structure  function  that  is  power  law  above  the 
inner  scale  <o.  The  exponent  in  (23)  is  large  unless  two  of  the  first  four  structure  functions  are 
small.  (The  other  two  can  nearly  cancel  the  last  two  structure  functions).  The  only  way  the 
exponent  can  be  small ,  while  allowing  the  cancellations,  is  for 

J'=0(so)  and  Jt  -  O(s0)  or  f'  =  O(s0)  (26) 


where  »o  **  the  field  correlation  distance  defined  by  Z?e(so)=  1-  The  remaining  variable  0'  or 
f )  is  typically  of  order  of  the  scattering  disk 


r  =  0^  = 


R 

*o* 


(27) 


s 


where  80  is  the  width  of  the  anguter  spectrum.  When  the  inner  scale  is  larger  than  W,  there  is 
little  scintillation.  Therefore,  we  consider  the  case  s0«CW,  y-  G(s0), and  £*  =  0(W).  A 

Taylor  series  expansion  about  ft  reduces  Q  to 

Q  =  +  (28) 

z 

For  a  power-law  structure  function 

=  (29) 

Q  becomes 

Q  .  E&L liL  ( =  (30) 

On  the  other  hand,  the  other  terms  in  the  exponent  are  O(D$(s0 ))  =  O(l),  which  was  the  condi¬ 
tion  which  caused  f  and  3’  to  be  order  So-  Thus,  if  Sq<W,  a  Taylor  series  in  Q  is  appropriate, 
but  a  Taylor  series  in  the  entire  exponent  requires  many  more  terms. 

The  other  possibility,  =  O(s0),  and  f-0(W)  requires  an  expansion  in  another  variable, 
O',  obtained  from  Q  by  interchanging  with  f.  This  alternate  expansion  is  the  fundamental 
reason  that  two  different  series  are  required. 

We  now  return  to  the  Taylor  series  expansion  in  Q,  with  the  result 

G<(  *;t\,7?)  =  £  Gf  (*;*../?)expf  ~  ±1  D9(f  +  S)  *  D9(f  -3»]  (31) 

nsO  2 

1  +  t  S-/  «.*•■/  **»exp  fW*/ Xcos^-.A;) -cosfr.*,)] 

-•  —  i  Ji=» 

. 

This  series  should  converge  quickly  when  the  quantity  Q  is  small  over  the  important  region  of 
integration.  Note  that  the  symmetry  of  the  full  fourth  moment  expression  (23)  does  not  hold 
for  each  term  of  the  expansion.  The  equivalent  moment-equation  derivation  is  presented  in 
section  4.1. 

The  first  term  of  (31)  reduces  to 


=  G!(tI.*e.rl.r8./f)G8(*8.  vr9,r4.^)  (32) 
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Thus,  (or  spatially  coherent  sources,  the  first  term  of  the  fourth  moment  is  a  product  of  two 
second  moments,  i.e. 

r4o(*l.ta.t3l*4./?) = )  (33) 

The  next  term  of  the  Green’s  function  expansion  is 

=  ZG{  (*;*<.*  )exp[  -  +  3’)  +  Z?*(f  -  3’)]]  (34) 

m 

/*e(<)[cos(3*.«)  -cos(f  •  *)]e‘^**  dt 

These  two  terms  contain  the  useful  first-order  description  of  the  fourth  moment.  The  rate  of 
convergence  of  the  series  is  determined  by  the  higher  order  terms.  For  the  plane-wave  case, 
the  fourth-moment  expressions  generated  from  (33)  and  (34)  are  identical  to  the  asymptotic 
results  of  Zavorotnyi  [  1979a  Eq  (7)]. 

2.4  Intensity  Correlation 

There  are  few  measurements  of  the  full  fourth  moment  of  WPRM  [Gurvich  et  al,  1978, 
1979a].  However,  the  intensity  correlation,  a  special  case  of  the  fourth  moment,  is  commonly 
observed.  We  will  now  demonstrate  that  our  fourth-moment  series  (31)  generates  two  different 
expressions  for  the  intensity  correlation  (4),  one  valid  at  low  spatial  frequency,  the  other  valid 
at  high  spatial  frequency. 

The  low  frequency  version  of  the  intensity  correlation  is  obtained  from  the  fourth 
moment  (3)  by  setting  1|  =  tg  and  9s  =  9«  or  by  setting  3  =  0  and  f  =  0.  The  n=0  term  produces 

cV(a.?.*)=  “jr/r.(a.?-.r.«'.o)exp(  -a).  ♦  tf  -#•)•?]]  (35) 

exp(  -  ♦  I')  +  D,!f  -  »)))  if  it  if  dt' 

For  spatially  coherent  sources 

cy  (&j.R)  =  </(  ^-.R)><J(  ^-.R)>  =  </(!,.*  )x/(t 8.*  )>  (38) 


li.V :.v.- 
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or 

*yo» -n.^)  =  +  §-.j?>x/(A-|-.>f»  «-**•*  df  (37) 

which  normally  describes  the  smallest  spectral  scale  of  the  process.  Indeed,  for  plane-wave 
conditions  =  </>23(4)  since  the  average  intensity,  </>,  is  a  constant 

The  corresponding  terms  for  n=l  may  be  written  as 

<Y(a.M)=  ^^7r4(a\/§\?\3\o)exP  ~i|-[(a- a*)*  3' +  <£-?')•?*]  09) 

exp  +  3*)  +^e0'“3')]  ♦eOWcostf*.*) -cos{f  .*)]e<?‘*  da*  d£’  d?  d3*  d2 

and 

*{/(d  4/?)  =  ^^r/r *(&’#.?=  -4f-.3,,o)exP[-i|r(a-a,)-3,)eil*-«^  09) 

exp  +  *)  +  ^eOf-  -3')]  *eO)[cos(3,.4)-cos(f'.jJ))  dd'  d£’  d3*  d* 

This  term  describes  the  refractive  focussing  by  irregularities  of  the  size  of  the  scattering  disk 
(Rickett  et  al,  19B4).  For  the  case  of  plane  waves,  (39)  reduces  to  the  familiar  result 

*i ;0M.* )  =  4#e(a)exp[  -  D9(^~)]sm2[q2~]  (40) 

which  also  describes  the  low  q  behavior  of  ♦(/! The  Born,  or  weak-scattering  approxima¬ 
tion  is  obtained  by  ignoring  the  exponential  term  of  (40).  Similarly,  the  Born  approximation 
for  the  general  case  is  obtained  from  (39)  by  ignoring  the  last  exponential  term. 

The  high  frequency  version  of  the  intensity  correlation  is  obtained  from  the  fourth 
moment  (3)  by  setting  i|  =  1,  and  =  I3  or  by  setting  3  =  0  and  f  =  0.  The  n=0  term  produces 

<#(*.?.*)  =  {^yir*(iA,^'r,i‘'0)exp  -<|-K*-*,)-i,  +  (W)-/*,l 

expl-  +  3  )  +  D%{f  -  3')J  d&‘  df  dl  dtt 


(41) 
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For  spatially  coherent  sources 

cV  (*!,*.*) = (!».*.*)  (42) 

This  expression  is  the  high  frequency  approximation,  i.e.  the  intensity  correlation  is  the 
square  of  the  mutual  coherence  function,  and  may  also  be  derived  by  assuming  that  the  com¬ 
plex  electric  fields  are  zero  mean  Gaussian  random  variables. 

The  n=  1  term  reduces  to 

cy (*.?.*)=  r4(a\2'.r.*'.o)exp  -<|-[(a -<*•)•$'♦  (?-?’>•/&•]  (43) 

exp  -~[DB(f  +  i’)+D9(7  -3*)]|*e(ie)(cos(^.it)-cos(r-e)]e<?  ,<  d&'  dfr  df  di'  of* 


^{•/(d.4^)=^^r7r4(a'ii=-i)^.f.3>.o)exp[-i|-(a-a')-3<]e "  *  *  (44) 

exp  -  —  [Dt(f  +  3‘)  +  D*{f  -3')]  $e(*)[cos(3'»2)  -cos(f.*)]  d&'  df  dt'  di 
z 


Rickett  et  al  [1984]  argue  that  the  physical  mechanism  for  this  term  is  the  modulation  of  the 
small  scale  structure  (42)  by  the  large  scale  refractive  process  (38). 


We  have  shown  how  the  fourth-moment  series  generates  two  expressions  for  the  intensity 
correlation:  one  (3  =  0.^  =  0)  useful  at  low  spatial  frequencies,  and  the  other  (3  =  0,^  =  0)  useful 
at  high  spatial  frequencies.  The  region  of  validity  depends  on  the  statistics  of  the  phase 
fluctuations  and  the  initial  source  distribution.  The  case  of  plane  waves  incident  on  a  random 
phase  screen  with  a  power-law  spectrum  has  been  investigated  by  Gochelashvily  and  Shishov 
[1975].  Their  calculations  of  the  first  few  terms  of  the  intensity  spectrum  [Figure  l)  imply  the 
two  series  converge  quickly  when  the  quantity  Q  is  small  over  the  important  region  of  integra¬ 
tion.  The  rate  of  convergence  is  difficult  to  determine  a  priori;  the  contribution  to  the  inten¬ 
sity  spectrum  from  the  higher  terms  is  the  best  indication  of  convergence.  A  finite  number  of 
terms  from  the  two  series  can  be  merged  by  a  weighted  sum  based  on  this  rate  of 
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convergence.  In  strong  scattering  conditions,  these  two  series  merge  quickly  and  only  a  few 
terms  are  required  to  describe  the  complete  spectrum.  The  intensity  correlation  series  is 
more  difficult  to  interpret  since  the  errors  of  the  expansion  accumulate  in  the  region  of  small 
spacing  but  in  the  spectral  domain,  these  errors  appear  in  the  central  regions  of  the  spec¬ 
trum  [cf  Figure  lj. 

3.  PATH  INTEGRAL  TECHNIQUES  FOR  EXTENDED  RANDOM  MEDIA 
3.1  Introduction 

We  now  consider  the  more  complex  problem  of  wave  propagation  in  a  random  media  that 
is  locally  homogeneous  with  statistics  that  vary  slowly  in  the  direction  of  propagation.  Laser 
propagation  in  the  atmosphere,  radio  propagation  through  the  interstellar  medium,  and 
sound  propagation  through  the  ocean  with  no  deterministic  background  are  common  exam¬ 
ples  of  this  phenomena.  Many  theoretical  methods  have  been  applied  to  these  problems. 
Moment-equation  methods  are  reviewed  by  Prokhorov  et  al  [1975].  Recently,  Macaskill  [1983] 
and  Frankenthal  et  al  [1984]  have  applied  the  two-scale  embedding  procedure  (Frankenthal  et 
al  [1982],  Beran  et  al  [1982])  to  produce  a  solution  for  the  fourth  moment.  The  application  of 
path  integral  techniques  to  problems  of  wave  propagation  in  random  media  was  introduced  by 
Klyatskin  and  Tatarskii  [1970],  Zavorotnyi  et  al  [1977]  and  Dashen  [1979].  Path  integral 
methods  have  been  successfully  applied  to  the  difficult  problem  of  WPRM  for  anisotropic,  inho¬ 
mogeneous  medium  with  deterministic  background  of  refractive  index  [Flattfe  et  al,  1979].  A 
functional  operator  form  of  the  path  integral  was  used  by  Tatarskii  and  Zavorotnyi  [i960],  to 
extend  thin-screen  analysis  to  the  problem  of  wave  propagation  in  extended  random  media 
for  the  plane-wave  case.  We  use  the  path  integral  representation  of  the  Green’s  function  to 
illustrate  this  connection  because  one  obtains  a  clear  presentation  of  the  role  of  the  source 
distribution.  We  believe  that  the  operator  formalism  is  equivalent  to  the  path  integral  method" 
and  produces  the  same  results. 

The  path  integral  technique  is  introduced  by  a  review  of  the  second-moment  derivation. 
The  reduction  of  path  integrals  to  familiar  Reimann  integrals  is  performed  by  a  useful  identity 
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[cf  (62)].  Using  this  identity,  we  present  a  series  expression  for  the  fourth  moment  that  is 
analogous  to  the  thin-screen  results  of  section  2.2.  The  behavior  of  the  resulting  correlation 
series  is  then  discussed. 

For  narrow  angular  scattering,  the  scalar  field  satisfies  the  parabolic  wave  equation 

ZUc^-  +  VzE  +  2JciTi(t,z)E  =  0  (45) 


where  n(i,z)  denotes  the  refractive  index  fluctuations.  We  define  the  correlation  of  refractive 
index  fluctuations.  Bn(&,t,z),  as 

5n(3.f.z)  =  <n.(0.z)n(a,z +f  )>  (46) 


The  path  integral  formulation  for  the  Green’s  function  was  developed  by  Feynman  [1948],  and 
may  be  written  as 


C(fAR)  =  f  Z7f(z)exp 


i§-/[^(z)]2dz  -i*/n[?(z).z]dz 
c  0  0 


(47) 


where  Dt{z)  denotes  the  infinite  dimensional  integration  over  all  possible  paths,  f(z),  con¬ 
necting  the  points  (f,0)  and  (i./?)  and  ?(z)  =  The  most  important  paths  are  those  near 
the  geometrical  path  from  (t*,0)  to  (t,/?),  given  by 


fc(z)=i'(l-”)+l|- 


(48) 


Transformations  of  these  geometrical  paths  will  be  denoted  with  the  subscript  G. 

If  there  are  no  refractive  index  fluctuations,  the  Green’s  function  becomes  the  free  space 


Green’s  function. 


S' (**./?)  =  f  £>f<z)exp 


c  o 


ZniR 


exp 


ik 
2  R 


(*-*)! 


(49) 


3.2  Second  Moment  by  Path  Integral  Techniques 

The  Green’s  function  for  the  moments  of  the  field  are  easily  expressed  in  terms  of  the 
path  integral.  The  Green’s  function  for  the  first  and  second  moment  were  derived  by  Dashen 
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[1979].  These  results  were  obtained  under  the  Markov  approximation  [Zavorotnyi,  1978].  We 
review  the  derivation  for  the  Green's  function  for  the  second  moment,  which  is  given  by 

G2(t„te;r1.W2./?)  =  <£(*,;*„/?)<?  (*e;*2./?)>  (50) 


=  //Z?fj(*)r>f2(*)exp  i|/([^(*)3*-[f2(2)]2)d*  <exp  -i*f[n(tl(z),t)-n(f2(z),z)]dz> 


Applying  (19)  and  the  Markov  approximation  results  in 


G2(ii;l'i,/?)=  Dt,(z)  Df2(z) 

I  u  R 

i  .•*  f/rX  m2  rl  i 


exp  if -f([^(z)]Z-[h^)]2)dz -^fd[tM-h(z).z]dz 
<-  *0  *  0 


and  the  phase  structure  function  density  is  given  by 


dtf,z)  =  &zf[Bn(0.t.z)-Bn(fi,t,z)]dt 


=  4rrfc2/ [  1  -  cos(fi  •  $)]$„  (4-9*  =  0»z  )  4 


where 


*n(4<f«./?)  =  ^^77fl>(t».*)exp[  ~<(t*$  +  z7«)l  dl  *  (53) 

is  the  refractive  index  spectrum.  Change  path  variables  to  the  centroid  and  difference  coordi¬ 


nates 


A(*)  =  7[*»(*)  +  *2(*)] 


^(z)  =  f,(*)-f2(*) 


The  double  path  integral  can  be  evaluated  by  expressing  the  paths  as  deviations  from  the  two 


geometrical  paths  defined  by  (48).  i.e. 


Ai(*)=A(*)-Ac(*) 

0i(*)=A(*)-0c(*) 


■  •  ■  >  ■  »  •  -  - 
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where 


Ac(*)sA(0)(l  -  +  *,)(  1  -  ^-)  ♦  |(t,  + 1,) i- 


Ac(* )  =/J(0)<l  -  %)+p(R)  “-=  (**,  -f8)(l  -  |-)  +  (i,  -  i*)  J- 
Then  Ai(°)  =  Ai(^  )  =  £i(0)  ~Pi(R)  ~  0  and  the  second-moment  Green’s  function  becomes 

GtOiXi.R)  =  exp  -l,)2  -  (*2  -  ta)*]  //  Z?Ai(* )  DPiiz ) 


>[*/A,(*M, 


(z)dz  ~-r/d[Ai(z)  +  £c(z  ).*]<& 
z  0 


Integrating  the  first  path  term  in  the  exponential  by  parts  and  substituting  the  free  space 
Green's  function  results  in 


R  R 

exp  -ikf  jii(z)>'fix(z)dz  -^fdifi^z)  +flc(z).z)dz 
0  *  0 


This  path  integral  is  evaluated  with  the  identity 

f  f  Dfi(z)  Dp(z)F[fi(z)]ex p  */>(*)•[?(*) -6(*)3«**]  =  (62) 


where  g(z)  is  the  solution  of  *£(*)  =  £(z)  that  satisfies  the  appropriate  boundary  conditions. 
Identity  (62)  reduces  (61)  to 


ChiK* t.ff)  =  O' (i,;*,./?)^  (***2.*  )exp  - dtfc(z ).z ]dz 

z  0 


and  J(z  )  =  0  is  the  solution  of  |(z  )  =  0  ,  that  satisfies  the  required  boundary  conditions. 


la 


3.3  Fourth  Moment  by  Path  Integral  Techniques 

Using  the  path'integral  technique,  we  derive  a  series  expression  for  the  fourth  moment  in 
a  manner  analogous  to  the  thin*screen  derivation.  The  Green's  function  for  the  fourth 
moment  is  given  by  the  multiple  path  integral 

<?«(*;*«.*)  =  (*4.*'4./?)>  =  /  Dtx  Dtz  Dtz  DtK 

R  R 

exp  » |*/[*?- >*  + *3  “fc «]d*  <exp[-ar/fn{Plt*)-n(f2,8)  +  n(fs.*)-n(f4.*)]cfa]>  (64) 

.0  0 

R 

«/[  n(t\.z)  -n(t 2.*)  +  n(?3,z)  -n(tA,z)]dz  is  a  zero  mean  Gaussian  random  variable  (19)  and 
o 

the  Markov  approximation  is  valid,  the  last  term  of  (64)  becomes 
R 

exp(  - 7- f[d(t  1  -  f2.z)  +  d(f3 - U.z )  +  d(?i  - ?4.z )  +  d(fz - f3,z )  - d(f2 -  ?4.z )  - d(t ,  - f3.z )]dz  ]. 

2  0 

Change  path  variables  to 

za(z  ) = ?,(z )  +-  ?2(z )  *-  /3(z )  ♦  ?4(* )  2  t,(z)=a(z)+fl(z)  +  j(z)+6(z) 

2 fi(z  )  =  fx(z  )  +  f2(z  )  -  f3<*  )  -  *4(*  )  2?2(*  )  =  a(z  )  +  0(z  )  -  ? (z  )  -  J(z  ) 

2f  (z )  =  t  »(z )  -  ?2(z  )  -  fs(z )  +  f4(z  )  2f,(z  )  =  3(z)  -0(z  )  -  f  (z  )  +  J(z  ) 

2i(z)  =  f»(z)-f2(z)  +  f3(*)-f4(z)  2f4(z)  =  a(z)-^(z)  +  f(z)-3(z)  (65) 

Following  the  second  moment  derivation,  we  express  the  paths  as  deviations  from  the 
transformed  geometrical  paths,  i.e. 

a,(z)=d(z)-ac(z)  (66) 

&  (*)=0(*)“0e(*) 
fi(*  )==?(*  )-tc(*) 

0i(*)=0(*)-0c(*) 

where 

ac=a(o)(i  -|-)+a(/?)f-  (6?) 
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In  order  to  simplify  the  large  expressions,  *]  will  replace  [d(0),£(0),f(0),3(0)]  and 

[a.#.?.3]  will  replace  (a(fl),jJ(i?),}(/?).3(/?)]  .  Then, 

GtiWuR)*  tejp-GUWi.R)/  DA,  Dfa  Df,  Z?3,exp  */(*, X  ♦Ml)**  (68) 


exp|-  ~.Ad(?i  +  3|  ♦  Jc(*  )+3c(*  ).* ) +  *(?  i  “  3,  +  ?c(*  )~3c(*  ).*  )  ♦  d(fa  +  3,  +  £c(e  )+3  c(z  ).*) 

+  d(fa  ~  fa  ♦  ?c(z  )-3c(2  ).* )  -  *  (&  +  f ,  ♦  &?(* )+?  c(*  ).*)-*  (fa  ~  ti  +  £c(*  )~f  c(*  ).*  )]**] 

Substituting  (integration  by  parts) 

R  R 

*/(4 r&i+tfh)*  ~  -Dcf(Arl,  +  fa>'i,)dz  (69) 

0  0 

and  performing  the  integration  over  the  paths  A,  and  using  (62)  results  in 


Gt(%;*i,R)  =  (*;*4J?)//  Dfa  0f,expj  -tt|J,.?,d* 

a  , 

exp  -  jjfid(K  +  f  c(«  )+3c(*  ).*)  +  * (7i  ♦  ?c(*  )~3c(*  ).* )  ♦  e(d (£,  ♦  £c(*  )+3c(*  ).* )  (70) 

+  d(fa  +  $c(*)~3c(e),*)  -d(£,  +  ^c(*)+fc(*).*)  "*(£i  ~  ?»+&?(*  )”?<?(*  ).*)]}**] 

Note  that  (70)  is  similar  to  (83)  of  the  thin-screen  derivation.  The  parameter  e  has  been  added 
to  label  the  quantity 
K 

Q=-jj£[d{fa+fa(z)+tc(z),z)+d(fa+$c(*)-lc(*)'*) 

i +^c(*  )+fc(«  ),  *  ) -d(#, ^(s  )-^c(s ),» )  jd* 

*  - 

=  4n*,^’/*w(49,  sO.*)exp[i^.(^,  +  |c(«)]Kco*(4*3c(*))-co*W-(f,  +  fc(*)))]  d$d*  (71) 


as  the  analogous  expansion  quantity  to  the  thin  screen  derivation.  For  a  point  source,  this 
expression  is  small  in  the  contributing  regions  of  path  space  [Dashen,  1979].  Performing  the 
Taylor  series  in  Q  results  in 


G4(tr.tV./?)=  I  Dfa  flf, 


exp  “  7  »(* ) +  f  c(*  )  +  d(%(z )  + 1  c(*  Hc(*  ).*  )]<** 


exp[-ifc/?i(*).^j(*)d2]+  £  enf  dz ,  •  •  •  f  dznf  d^  •  •  f 

0  n=l  ‘  0  0  -«• 


exp  i(*)‘[iM*)-£  ~zi)]dz  +  i££c(z>)*fy 

,o  i- 1  * 


ft  (^j  -9*  _  0.z>  )[  COS  (^•3c(*i))  -cos 
1=1 


Applying  identity  (62)  yields 


=  exp  -~f[d(fc(z)+ic(z).z)  +  rf(fc(*)-^c(* ).*)]*  * 


£  *,.../  dZnf  rfq, . . .  /«/*, 

n*l  n!  0  0  —  — 


exp[t£flc(zj)*%  "“■/[<*(  £  %-fc(z,.*m)  +  M*)+$c(* ).*)  +  «*(  £  %-h(2;xm)  +  fc(z)-^c(z).z)]^ 

l  j  =  l  z  0  msj  *  msi  *  , 

ft*»(^.9#  =  0.«j)[cos(^.Jc(*f»  ~cos(^ .(  £  +  ^c(*i»)]l  =  £  CAn(ii:ti,R)en 


where  A(z;Zi)  is  the  solution  of 


A(z;*,)  =  d(*  -*t) 


that  is 


A(s;s,)  =  *(-^--l)  *<*, 


h(e;*,)  =  *,(^--l)  *>*, 


v‘*  -V  *m  *«  _  *«"  •«*  *  *  •_*  •'  .  ’  .  •  .  * 
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The  free  space  Green's  function  for  the  fourth  moment,  G(  (*(;?<.#).  is  given  by  (24).  The 
terms  of  the  series  are  identified  as  the  coefficients  of  the  parameter  c*.  This  construction 
will  be  used  to  prove  the  equivalence  of  the  path-integral  series  to  the  iterated  moment- 
equation  series  of  section  4.2 
The  first  term  of  (73) 

M 

-  ♦*<?<•  >•»>  ♦<»«(•)-*«(*  >•*)!*• 

C4o(Ii:fi.f?)  =  C{(tt;4'i./?)e  0  (76) 

is  analogous  to  the  thin-screen  result  (32). 

The  n=  1  term  is  more  complex. 

*  •  , 

G.(ii;*i.R)  =  4nkzG{(1ti-.ri.R)f f4n($<l,=(>'zi)[c°s(A-&c(zi))-cosfa-[f~h(zl.zl)  +  yc(zt)])] 

o  —  K 

e‘4-?c(*,)exp  _i_  f[d(^h(z.z  ,)+?<;{*  )+M* ).*  )+<* ( I’M* ;* I )+?c(*  )~M*  ).* )]dzj<i^dz  1  (77) 

The  higher  order  terms  of  the  fourth-moment  Green's  function  are  obtained  from  (73) 
but  become  more  intractable.  However,  the  first  order  description  of  the  fourth  moment  is 
given  by  the  leading  terms  of  the  series  expansion.  We  now  compare  these  expressions  to  pre¬ 
vious  results  for  the  fourth  moment.  The  point-source  result  is  obtained  from  the  Green's 
function  by  setting  to  zero.  This  reproduces  the  iterated  series  of  Shishov  [1972  Eq  (13)] 
that  was  derived  using  moment-equation  methods  in  a  spherical  coordinate  system.  Appling 
the  plane-wave  initial  condition  generates  the  series  expression  of  Shishov  [1971  Eq  (20)]. 
Gurvich  et  al  [1979b]  proved  this  series  was  convergent.  It  can  be  shown  that  the  point- 
source  series  is  also  convergent.  Applying  the  plane-wave  case  to  (77)  with  d  =  0  and  2  =  0, 
reproduces  the  strong  scattering  results  of  Fante  [1975  Eq  (Bl)]  and  Zavorotnyi  [1979b  Eq 

(3)). 
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3.4  Intensity  Correlation 

Furutsu  [1972]  derived  expressions  for  the  intensity  correlation  from  a  Gaussian  beam 
propagating  in  a  random  medium  with  a  square-law  structure  function.  This  case  describes 
random  wander  of  a  beam  [Wandzura,  I960].  Weak  scattering  results  for  intensity  correlation 
from  a  Gaussian  beam  are  presented  by  Ishimaru  [1969].  We  now  consider  the  intensity 
correlation  following  the  thin  screen  analysis  of  section  2.3. 

The  fourth-moment  series  generates  two  expressions  for  the  intensity  correlation;  one  is 
obtained  by  setting  3=0  and  f  =  0  (the  low  spatial  frequency  region)  and  the  other  by  setting 
3  =  0  and  $  =  0  (the  high  spatial  frequency  region).  The  resulting  expressions  for  the  n=0  term 
are 

4!  -*)•* (?8) 

exp  +  3*)(l--|-).*) +  <*((?' -*')(!- ~).z)]dz  d&  dpdfdi' 

and 


r4(a\*\t.3\0)exp[  -i|-[(a  -  *•)•*  +  (?  -?).*-)]  (79) 

exp[-^-/[d((r  +  3*)(l-|-)  +  7|-.*)  +  ‘i((r-i,)(l-|r>+f|-.*)]^  dvdpdfdt' 

For  spatially  oherent  sources,  these  expressions  reduce  to 

cy(aj,tf)  =  </(^-./?)></(~^,/?)>  =  </(*„  2)></(^.2»  (80) 


and 

<#(*,.*2,/?)  =  r2(*,.*2./?)ra*  (*,.*2./?)  (si) 

which  are  analogous  to  the  thin-screen  results  (36)  and  (42). 

The  corresponding  expressions  for  the  n=l  term  are  given  by 


expj-  ±-|[d(  ;t ,)  +  (f  ♦  3*)<1  “  jfU ) ♦  i)  ♦  (r  -3‘)(  1  “  f-).*  )}«** 


[cos(^3*(i  -  ^)-cos(4.(  -  ^-))»  df  d*  dA**x  (02) 


c^/(d.f./?)=^^|/r4(a^r.3,.°>e  *  c 


*  *n(**l) 


jcos[«l-3'(l-~^]-cos[il.(J-A(zl,z,)  +  f(l-^-)  +  f^-))  dV  dp  df  dl  dUdzt  (83) 

The  appropriate  Fourier  transform  (cf  (5))  produces  the  corresponding  expressions  for  the 
intensity  spectrum.  The  Born  approximation  is  secured  from  (82)  by  ignoring  the  last 
exponential  term. 

The  low-frequency  series  for  the  intensity  spectrum  has  been  investigated  for  a  power-law 
spectrum  of  refractive  index  fluctuations  and  plane-wave  [Gochelashvily  and  Shishov,  1971, 
1974]  and  point-source  [Gochelashvily  et  al,  1974]  geometries.  The  qualitative  behavior  of  the 
intensity  spectrum  is  similar  to  the  thin  screen  result  [Figure  1].  In  weak  scattering,  the  first 
two  terms  describe  the  complete  intensity  spectrum.  In  strong  scattering  the  intensity  spec¬ 
trum  is  characterized  by  two  components.  The  low  frequency  region,  ♦/(jj./?),  is  due  to 
refractive  focussing  by  irregularities  of  the  size  of  the  scattering  disk  (Rickett  et  al.  1984). 
The  high  frequency  behavior  was  first  determined  by  summing  the  low-frequency  series,  which 
yields  the  high  frequency  approximation  ♦o The  derivation  presented  here  generates 
this  same  result  plus  higher  correction  terms  without  performing  a  complicated  summation. 
As  in  the  thin  screen  case,  the  physical  mechanism  for  the  first  correction  term  (83)  is  the 
modulation  of  the  small  scale  structure  (81)  by  the  large  scale  refractive  process  (82). 


Uscinski  [1982]  derived  an  approximate  expression  for  the  intensity  spectrum  for  the 
case  of  plane  waves  by  summing  a  perturbation  series  for  the  fourth  moment  equation. 
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Macaskill  [1963]  derived  the  same  result  using  the  two-scale  expansion  [Beran  et  al.  1962].  In 
our  notation,  their  expression  for  the  intensity  spectrum  is 

*(4*)=  (^jr/exP  <*(*-*)+  d(f-(* -*).«) 

e-1'**  d*  (84) 

The  extreme  low  and  high  spatial  frequency  behavior  is  identical  to  that  obtained  from  (82) 
and  (81).  However,  the  predictions  for  the  intermediate  frequency  region  are  different.  More 
theoretical  calculations  are  required  in  order  to  determine  the  accuracy  of  the  two  methods. 

The  intensity  correlation  from  a  point  source  embedded  in  a  random  medium  with  irre¬ 
gularities  that  are  constant  in  one  transverse  dimension  is  obtained  from  the  point-source 
result  by  the  substitution  $n(«D  =  $n(?*)^ (9y)>  The  leading  order  terms  of  the  intensity  corre¬ 
lation  series  are 

=  </ >*exp[  )  (85) 

/?  «•  2 

CY  (P.R)  =  **k*<I>2ff*n(q.zl)s\nz[lg-h(zuzl)] 

X  z 

exp  -/d(  j-ft(z;zi),z]dz  +  ipq~  dq  dzt  (86) 

and 

CVifi.R)  =  Qnk2<I>*ffin(q  ,z  ,)sinz[g  p  +  ^-h(z„z ,)] 

x 

exp  -jTd[Pjj-  *  J-Mz;«,).*]dz  +iPq~jf  dq  dzt  (87) 

Ignoring  the  exponential  term  of  (86)  reproduces  the  weak  scattering  (Rytov  approximation) 
result  of  Tur  and  Beran  [  1983]. 


Fante  [  1983]  investigated  the  effect  of  the  inner  scale  of  turbulence  on  scintillation  for 
the  case  of  strong  scattering  of  plane  waves  by  calculating  the  expression  C*/ (£./?)•  The 
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contribution  of  this  term  to  the  total  intensity  variance  is  appreciable  over  a  wide  range  of 
parameter  space.  Therefore,  higher  terms  are  required  for  an  accurate  description  of  the 
intensity  correlation. 

4.  MOMENT  EQUATION  APPROACH 

4. 1  Moment  Equation  Method  for  the  Thin  Screen 

We  now  derive  the  fourth-moment  result  (23).  using  the  moment-equation  method.  This 
method  is  based  on  partial  differential  equations  for  the  transverse  moments  that  are  derived 
from  the  parabolic  equation  for  the  random  fields.  These  differential  equations  may  then  be 
solved  by  transform  methods  [Rumsey,  1975]  and  the  method  of  characteristics  [Kiang  and 
Liu,  1982],  There  are,  in  general,  many  transforms  that  will  simplify  the  problem.  We  choose 
one  [Shishov,  1971]  that  permits  an  analogous  derivation  for  the  fourth  moment  of  waves  pro¬ 
pagating  through  extended  random  media. 

The  fourth  moment  satisfies  the  differential  equation 

^-(V?-V|  +  V|-V|)r4(*1.Je.t3.t*.*)  (88) 

with  initial  condition 

r4(*i.*2.*3.*».o + )  =  = r4(t,.tz.t3.»4.o) 

exp[  -  ±-[ 0*(i,  -  **)  +  z>e( la  -  *.)  +  /?«(*,  -  *)  +  Dfa  -  *,)  -  D6(h  -  *,)  -  />«(*,  -  *g)]](89) 
Change  variables  to  the  coordinate  system  (22).  Then 

^-(M,*3.*)=  -(V  v, + Vj.  Vj)r4(a.j&.f.$.* )  (90) 

This  equation  is  solved  in  the  Fourier  transform  domain 

Af(*,4f.a./?)=  ^jr/n,(aAf.3.j?)e-<i*-«^>ida  dfi  (91) 


which  transforms  (90)  to 
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fjf (92) 

Changing  variables  to 

&  =  jfc8 -<z  6=  fc}  -  $z  r  =  z 

J(r)=irr- 

V,  =  t»,  V*'»  (93) 

transforms  (92)  to 

^■(*.aa.6.r)  =  0  (94) 

Since  M  is  independent  of  r,  the  solution  in  the  original  variables  is 

M  =  M(*Af  - -.3  - 2  —.0  +  )  (95) 

The  fourth  moment  is  the  inverse  Fourier  transform  of  this  expression,  i.e. 

rt{a.M.R)  =  7 —zrj  k  ,o)ei[i '(*-*">  *  *’<*'*'»  (96) 

(2^) 

exp  - ~[f?e(r  +  3’)  +  09(?  - 3')  +  D9$'  +  3')  +  Z?e(jfr  - 3')  - Z>9(/}’  +  f)  -  D9@  -?))  d*  dfi'  d^di 

The  change  of  variables  k  and  3’=3-jf/?/fc  produces  the  Green’s  function  result. 

(23).  The  moment-  equation  method  is  based  on  the  formulation  of  the  free  space  Green’s 
function  in  the  Fourier  transform  domain.  This  transform  technique  will  now  be  applied  to  the 
more  difficult  problem  of  wave  propagation  in  an  extended  random  medium. 

4.2  Moment  Equation  Method  for  Extended  Random  Media 

The  moment-equation  method  was  a  major  contribution  to  the  theory  of  wave  propaga¬ 
tion  in  extended  random  media  (Prokhorov  et  al  1975,  Tatarskii,  1971).  Urng  moment- 
equation  methods,  a  series  expression  for  the  fourth  moment  was  derived  for  plane-wave 
[Shishov,  1971)  and  point-source  [Shishov,  1972)  conditions.  For  an  arbitrary  source 
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distribution,  we  present  an  analogous  expression  that  is  identical  to  the  path  integral  results 
of  section  3.3,  thus  demonstrating  the  equivalence  of  the  two  methods. 

The  transverse  moments  of  the  electric  field  satisfy  differential  equations  [de  Wolf,  1967; 
Brown,  1967;  Shishov,  1968;  Dolin,  1968;  Chernov,  1969;  Tatarskii,  1969;  Lee,  1974].  We  follow 
the  techniques  of  Shishov  [1971]  to  solve  the  differential  equation  for  the  fourth  moment 

^[*?-vf+v|-vUr4+Kr4=o  (9?) 

where 

2V  =  d(il-tz.*)  +  <l(*i-*4.z)  +  d(ia-ta,z)  +  d(%i-i4,z)-d(ty-i2,z)~d(1t2-i4.z)  (98) 

In  the  coordinate  system  of  (22),  the  fourth  moment  satisfies 

f~r  -  f  [ ' • Vt  +  ] : r4 +  v  r4 =  o  (99 ) 

with 

2  K  =  d(?  +  3,f?)  +  d(?-3,./?)  +  e[d(£  +  3./?)  +  d(j&-3./?)-d(£  +  ^/?)-d(#-?.f?)]  (100) 

where  the  parameter  t  identifies 

5  =  ^jd(/5  +  «.2)  +  tf(^-3.*)-rf(^  +  f.e)-«i(?-f.2)]  (101) 

as  the  same  combination  of  structure-function  densities  we  used  in  the  path-integral  expan¬ 
sion  variable  Q  of  (71).  The  Fourier  transform  (91)  converts  (99)  to  the  integral  equation 

+  f  •  +  f  •  V*  +  (?  +  3)  +  *(?  -  Z))M  =  )  (102) 

where  G  is  the  convolution  of  M  with  the  transform  of  S,  Le. 

«a 

G(2,^,? ,3,2 )  =  4tt*2/*„(4j,2)[cos(3 . $i)  -cos(f  3i.?.3.z )d$i  (103) 

The  change  of  variables  (93)  transforms  (102)  to 

~  ♦  +  3 (r).r)  +  d(f(r)  -3(r),r)]A/  =  sGMf(r)J(r).r) 


(104) 


which  is  analogous  to  (94)  of  the  thin  screen  derivation.  This  is  an  ordinary  differential 
equation  in  r,  with  solution 


W(i.4b,S.r )  =  f,0)exp  -  ±f[d (f (2 )  +  3(z  ).* )  ♦  d  (f (z )  -  3(z  ).z  )]dz 


+  e/exp  -  —  j[d(? (z  )  +  3(z ).z )  +  d(?(z )  -  3(z ),z )]dz  G(*.4? (*i).3(z ,).* ,)  dx ,  (105) 


Changing  variables  back  to  3 ,f,R  results  in 

[  *  ) 

M(2,4?,3./?)  =  A/(*.4?3.33.0)exp  /[d(f2  +  32.z)  +  d(?2-32.z)]dz 


where 


+  £jfexp  -7'jT[d(f2  +  32lz)  +  d(^2-32.z)]dz  G(«.4fi.3,.zi)  dz,  (106) 


f,  =  f-|-(/?-zl)  fz  =  f-|-(/?-*)  ?3  =  ?~f/? 


3.  =  3- f-(tf-*i)  32=3-|-(/?-z)  33  =  3-|-/? 


(107) 


This  expression  can  be  written  as  the  integral  equation 

M(t.W,S,R)  =  Z(il.M.lR) 


X  -  -  - 


+  e/// /K(*,44\f,?\3,3\/?.z i)tf(X.4.f .3*.z ,)  d*  df  d3#  dz , 


(108) 


where 


^(*.4^.^.^)  =  ^(^.4^3.23.0)exp  --^-/[d(?2+32,z)  +  d(f2-32.z)]dz  (109) 

*  0 


9 


K(*.44.f  .f  .*.3'./?.*  i) =4»r***n  (4- *.*.)*(*'  -3,)3(f -1 ,) 
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(no) 


exp 


-  -“/[dto  +  *?•*)  +  dfc  "  *«•*  >ld* 


jcos[  J*  •  (4-  4)1  -  cos[r .  (4-4)]) 


The  formal  solution  is 


Af(jt,4,?,3,/?)  =  2  wn(e,4T, «,/?)£' 

n=0 


where 


(in) 


B  »  «  • 

Af„u.4M^)=// / / K(*.4«.f.r.*.*,.^.*i)Mi-.(«.4.t^,.*i)  <*$  *r  ** 


(112) 


and  4^.2, R)  =  Z(it, 4.?,3,/?).  This  solution  is  a  power  series  in  e.  The  path-integral  expan¬ 
sion  (73)  is  also  a  power  series  in  the  same  parameter  i.  Since  this  problem  has  a  unique  solu¬ 
tion  and  since  the  path-integral  and  moment  equation  are  equivalent  representations  (Codona 
et  al,  1985),  the  two  series  are  equal,  term  by  term.  We  will  show  this  equivalence  explicitly  for 
the  n=0  term.  The  fourth  moment  is  given  by  the  inverse  transform  of  (1 10),  i.e. 


exp 


-- jf[d(?+3-(4+*)(/?  -z )/  *  ,z  )+d(f -3-(4-*)(/?-e  )/  k  ,z  )]dz 


ddd/}'d4d* 


The  change  of  variables  ^-^-^R/  k  and  Z'-Z-itR/k  produces  the  same  result  obtained  from 
the  n=0  term,  (76),  of  the  path-integral  derivation.  Unfortunately,  for  higher  order  terms,  the 
equality  of  the  two  expansions  is  less  obvious  since  the  functional  form  of  the  two  series  is 
different:  the  moment  equation  series  is  essentially  a  multiple  convolution,  while  the  path 
integral  result  is  not.  However,  with  careful  algebraic  manipulation,  the  equality  of  the  two 
series  can  be  explicitly  demonstrated. 

The  moment-equation  method  was  first  applied  to  the  simple  geometries  of  plane  waves 
and  point  sources.  Since  the  moments  are  then  independent  of  the  centroid  coordinate,  the 
differential  equations  simplify.  Early  work  concentrated  on  moment-equation  methods 
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because  of  this  simplification.  Moment-equation  techniques  have  also  been  used  to  investigate 
the  validity  of  the  Markov  approximation  [Klyatskin,  1969;  Klyatskin  and  Tatarskii,  1971]  and 
the  effects  of  non-gaussian  refractive  index  fluctuations  under  the  Markov  approximation 
[Klyatskin,  1975).  The  Green's  function  formulation  provides  a  clear  connection  between  the 
thin  screen  and  extended  medium.  The  operator  form  of  the  path  integral  has  also  been  use¬ 
ful  for  evaluating  the  corrections  to  the  Markov  approximation  for  the  higher  intensity 
moments  [Zavorotnyi,  1978]. 

5.  SUMMARY 

A  series  expression  has  been  derived  for  the  fourth  moment  of  waves  incident  on  a  phase 
screen  or  propagating  through  extended  random  media.  These  results  can  be  derived  using 
moment-equation  techniques  or  functional  methods  (path  integral  or  operator).  The  asym¬ 
metric  terms  of  the  expansion  generate  two  expressions  for  the  intensity  correlation;  one 
that  approximates  the  low  frequency  region  of  the  spatial  spectrum  and  the  other  appropri¬ 
ate  for  the  high  frequencies.  The  rate  of  convergence  of  the  two  approximations  can  be  used 
to  produce  a  complete  expression  for  the  intensity  spectrum  valid  for  any  initial  source  dis¬ 
tribution.  The  calculations  required  for  a  complete  expression  may  be  excessive.  However,  in 
strong  scattering  conditions,  the  leading  order  behavior  of  the  intensity  spectrum  is  well 
described  by  a  few  terms  of  the  series.  The  expressions  presented  here  are  applicable  to 
many  problems  of  WPRM  that  involve  arbitrary  source  distributions.  These  include 

a)  The  intensity  statistics  from  beamed  lasers,  navigational  beacons,  radar,  spacecraft  and 
satellite  radio  transmissions,  astronomical  sources  and  other  extended  wave  sources. 

b)  The  effects  of  slowly  varying  refractive-index  statistics  along  the  propagation  path,  e.g.  the 
turbulence  profile  of  planetary  atmospheres. 

c)  Comparison  of  thin-screen  and  extended-media  results. 

d)  The  relationship  between  intensity  statistics  and  the  refractive-index  spectrum. 

e)  Imaging  through  random  media. 
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figure  Captions 


Figure  1.  The  leading  terms  of  the  intensity  spectrum  versus  normalized  spatial  fre- 

quency,  qRf,  where  /?/=(/?/*)*  is  the  Fresnel  scale.  The  curves  are  calculated 
from  expressions  given  by  Gochelashvily  and  Shishov  [1975]  for  the  case  of  plane 
waves  incident  on  a  random  phase  screen  with  a  Kolmogorov  spectrum  of  phase 
fluctuations  and  D%{Rj)= 100.  The  (-)  sign  indicates  that  ♦^(g)  is  negative  at 
high  frequency. 
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Abstract 

The  intensity  cross-spectrum  (spatial  Fourier  transform  of  the  two-frequency  intensity 
correlation)  for  scintillations  caused  by  a  plane  wave  passing  through  a  random  phase 
screen  is  considered.  Two  series  solutions  (one  valid  for  low  and  the  other  for  high  spatial 
frequencies)  are  obtained  which  are  the  generalizations  of  previous  results  for  the  mono¬ 
chromatic  intensity  spectrum.  We  show  that  the  Gaussian-field  approximation  (modelling 
the  cross-spectrum  as  the  transform  of  the  square  of  the  second  moment)  breaks  down 
when  the  outer  scale  is  large  compared  with  the  diameter  of  the  scattering  disk. 
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1.  Introduction 


Wave  propagation  through  random  media  gives  rise  to  intensity  fluctuations 
that  are  wavelength  dependent.  Examples  of  this  phenomenon  are  chromatic  stellar 
scintillation  [Jakeman  et  al,  1978],  pulsar  scintillation  [Rickett,1969j,  interplanetary 
scintillation  of  compact  radio  sources  [Cole  &  Slee,1980],  and  multifrequency  laser  pro¬ 
pagation  [Gurvich  et  al,  1979;  Azar  et  al,  1985].  In  weak  scattering,  the  intensity 
fluctuations  are  correlated  over  a  wide  range  of  frequencies.  However,  in  strong  scatter¬ 
ing,  the  intensity  fluctuations  are  decorrelated  after  a  relatively  small  change  in  fre¬ 
quency  [Popov  &  Soglasnov,  1985;  Cordes  et  al,  1985].  Further  complexity  arises  when 
the  intensity  correlation  between  spatially  separated  points  are  considered,  and 
different  wavelengths  at  these  separated  points  are  allowed  [Azar  et  al,  1985].  Finite- 
bandwidth,  finite-aperture  receivers  are  examples  of  such  cases.  These  problems  may  be 
analyzed  by  studying  the  two-frequency,  two-point,  intensity  cross-correlation,  or  its 
Fourier  transform:  the  two-frequency  intensity  cross-spectrum. 

A  common  approach  for  dealing  with  intensity  correlations  in  strong  scattering 
conditions  is  to  argue  that  the  real  and  imaginary  parts  of  the  field  are  the  sum  of 
many  independent  contributions  and  therefore  become  zero-mean  Gaussian  random 
processes.  We  call  this  the  Gaussian  field  (GF)  approximation.  The  GF  approximation 
implies  that  the  correlation  of  intensity  fluctuations  is  equal  to  the  magnitude  squared 
of  an  appropriate  second  moment  of  the  field.  In  particular,  the  two-frequency  two- 
point  intensity  cross-correlation  would  be  the  magnitude  squared  of  the  two-frequency 

second  moment,  which  contains  the  factor  exp[— -£-(A o/o)7#1],  where  a  is  the  center  fre- 

quency,  A  a  is  the  frequency  difference,  and  $  is  the  rms  phase  shift[Ginzburg  &  Erukhi- 
mov,  1971],  For  many  practical  experiments  $  is  large,  making  the  predicted  intensity 
decorrelation  bandwidth  extremely  small,  and  inconsistent  with  experiment  [Gurvich  et 
al,  1979;  Flatte',  1983].  Theoretical  approaches  in  the  1970’s  gave  various  reasons  for 
neglecting  the  rms-phase-shift  factor  completely  [Shishov,  1973;  Lee,  1976;  Dasben,  1979] 
Dashen’s  argument,  which  applies  to  both  the  thin-screen  and  extended-medium  cases, 
hinges  on  the  relative  size  of  the  outer  scale  (largest  sized  medium  fluctuations)  to  the 
diameter  of  the  scattering  disk.  He  argues  that  if  the  diameter  of  the  scattering  disk  is 
much  smaller  than  the  outer  scale,  then,  for  small  enough  frequency  differences,  the 
phase-shift  factor  can  be  ignored. 
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In  this  paper,  we  analyze  intensity  decorrelations,  in  frequency  and  transverse 
spatial  separation,  of  a  plane  wave  propagating  through  a  random  phase  screen.  A 
series  representation  of  the  two-frequency  intensity  cross-spectrum  (spatial  Fourier 
transform  of  two-frequency  intensity  correlation)  is  derived.  An  approximation  to  the 
cross-spectrum  is  presented  that  is  valid  for  all  values  of  the  outer  scale.  These  results 
extend  the  thin-screen  results  of  Dashen  to  arbitrary  frequency  differences  and  provide 
a  description  of  the  transition  from  a  small  (relative  to  the  outer  scale)  to  a  large 
scattering  disk.  The  r ms- phase-shift  factor  is  replaced  by  a  new,  more  accurate,  factor, 
and  the  conditions  under  which  this  factor  may  be  omitted  entirely  are  detailed.  In 
Section  2  we  introduce  the  second  and  fourth  moments  of  the  field  and  the  two-point, 
two-frequency  (intensity)  correlation  function.  In  Section  3,  we  consider  the  Fourier 
transform  of  this  correlation  function,  the  two-frequency  intensity  cross-spectrum.  We 
derive  two  different  series  representations  for  this  cross-spectrum.  The  first  series 
describes  the  low-spatial-frequency  behaviour  (i.e.  it  converges  quickly  in  this  region). 
The  other  series  describes  the  high-spatial-frequency  behaviour.  These  series  extend  (for 
the  thin  screen)  the  results  for  the  monochromatic  intensity  correlation  [Codona  et.  al., 
1985]  to  two  frequencies.  In  section  4  we  discuss  the  relationship  of  our  approximation 
to  the  GF  approximation.  It  is  shown  that  the  GF  approximation  (with  the  mean- 
square  phase  shift  factor)  is  valid  only  over  a  negligible  portion  of  the  intensity  cross- 
spectrum.  Finally,  in  Section  5  we  summarize  our  approximations  to  the  two-frequency 
intensity  cross-spectrum. 

2.  Definitions  and  Notation 

We  consider  plane  waves,  normally  incident  on  a  phase  screen,  that  freely  pro¬ 
pagate  a  distance  R  beyond  the  screen  to  an  observation  plane.  The  refractive  index 
fluctuations,  /i(x*),  in  the  screen  induce  a  random  phase  change,  0(x)=kf/j(jT,z)dz,  as 
the  field  passes  through  the  screen.  Here  z  is  the  direction  of  the  propagating  wave,  x*  is 
the  co-ordinate  transverse  to  this  direction,  and  k  is  the  wavenumber.  The  phase  change 
©  is  assumed  to  be  a  zero-mean  Gaussian  random  variable  with  homogeneous  statistics. 
We  consider  wave  propagation  which  is  characterized  by  narrow  angular  scattering  by 
the  phase  screen  due  to  the  small  fluctuations  of  the  refractive  index.  The  complex 
scalar  wave  field  can  be  expressed  as  E(jf,r ;*)«'**,  where  the  field  E  has  the  value  on 
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emerging  from  the  screen 


£(x*,0+;*)  =  £(x*,0;*)  e,0<r) 


expressed  in  terms  of  the  incident  field  £(jT,0;it).  In  the  space  after  the  screen  the  field 
satisfies  the  parabolic  equation 

dE  _  J_  &E_  (t>, 

dz  2k  dx*2  {  ' 

For  a  plane  wave,  the  incident  field  is  a  constant,  which  we  set  to  unity.  The  solution 
of  the  parabolic  equation  at  the  observation  plane  with  the  proper  initial  condition  is 


k  t  ,ikbC-x')2 


)  rfV 


Averages  of  the  fields  are  performed  by  the  use  of  the  identity 

—  —<©*> 

<«,e>=e  2  (4) 

valid  for  any  zero-mean  Gaussian  random  variable.  This,  for  example,  gives  the  average 
of  the  field  (first  moment)  as 


<£(x*  ,R-,k)>  =  exp  — 
£ 


where  &  is  the  mean  square  phase  shift. 

The  random  nature  of  the  fields  is  conveniently  described  by  statistical  moments 
of  the  field  evaluated  in  the  observation  plane.  The  moments  of  particular  concern  for 
the  study  of  frequency  decorrelation  are  the  second  and  fourth  moments.  The  general 
second  moment  is  given  by 

T^t[,r2,R\kuk2)  s  <£(x*1,/?;<r1)£*(3r2,/?;<r2)>  (6) 

^^2  ft f  /  ^  ■  i.  n  •  I  \  1 f  i  f  l  ^r*l\  2  I.  /-r*  J5_  I  ]2  I 


(2  xRf 


//exp  -(— -  kik2D„{rx '  -x*2' )  exp  ]rfVrfV 


where  the  fields  are  given  by  (3).  The  identity  (4)  is  used  to  obtain 


-J  <l/(  * iMno -*2M(r9rO  )*'|2>  __  (j)2*2 + * ,*aD<1{rl-ra) 

C  a  t> 


with  the  wavenumbers  expressed  in  terms  of  sum  and  difference  variables 
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f  =  (8) 

The  second  moment  is  expressed  in  terms  of  <J>a,  which  is  the  mean-square  phase  shift  for 
a  wave  with  the  mean  wave  number,  IF,  and  the  phase  structure  function,  k^D,,,  asso¬ 
ciated  with  a  wavenumber  which  is  the  geometric  mean.  The  wavenumber  dependence 
of  the  phase  structure  function  is  explicily  displayed  by  expressing  it  in  terms  of  a  struc¬ 
ture  function  for  the  integrated  index  of  refraction  fluctuations, 

Zyx-j-x-j)  m  <  d*T  >  (9) 

It  is  convenient  to  express  the  second  moment  in  terms  of  sum  and  difference  variables 
for  the  co-ordinates 

5*=  *‘*X-  ,  0  =x*i-x*2  (10) 

with  similar  expressions  for  the  (primed)  co-ordinates  on  the  screen.  Then  the  second 
moment  is  given  by 


Sk  VO  1  f  t  ok  >*9®*$"*) 


exp  MP{f  (11) 

The  a1  integration  is  immediate  since  the  phase  r2  is  independent  of  which  follows 
from  the  translation  invariance  of  the  problem  for  an  incident  plane  wave.  Performing 
the  a'  integration  we  secure  for  the  second  moment 


„  ,■?  „  .  .  »  *i*a  r  ,»***,  r  ~i kik*DAf*)  ik\k%  -► 

ra(/i,fi;M1)-^«p(-(T)1Tl/«  <12> 


The  general  fourth  moment  involves  the  field  at  four  different  spatial  points. 
From  the  solution  for  the  field,  (3),  this  moment  also  involves  the  integration  over  four 
points  on  the  screen.  It  is  convenient  to  transform  the  co-ordinates  of  the  screen  to 


1111* 

!  1  1-1-1  *2' 

7  1-1-1  1  rs' 

1-1  1-1 
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A  similar  transformation  exists  for  the  co-ordinates  of  the  observation  plane.  For  the 
two-frequency,  two-point  intensity  correlation  function,  which  is  a  special  case  of  the 
fourth  moment, 

C/(* s  </(x'l,J?;*l)/(xyf;*9)>  (14) 

~<E(?l,R\kl)Et(Tl,R-,kl)E(ri,R\k2)E'(?2,R;k2)  > 


two  of  the  transformed  quantities  on  the  observation  plane,  7  and  6,  are  zero.  Using 
the  solutions  (3),  the  above  transformation,  and  (8),  we  find  that  the  integration  over 
the  centroid,  a',  gives 


2*  m 


_±v 

<iv* 

[  e  *  exp 


(15) 


where  P  =x^  -x^,  and 

V4  =  Jfc?  D,(  (1  h  ' )  +  *|  D„(  (1  )7 ') 

2*  2k 

—k1k2\DM{P'+'f  ^DJP'-I  7  01 


(16) 


3.  Intensity  Cross-Spectrum 

The  intensity  cross  spectrum  is  given  by  the  p  Fourier  transform  of  (15),  which 

produces  the  delta  function,  —  7f(  1  —{6k/2t)2)),  where  T  is  the  transform  vari- 

H 

able.  Integration  over  7'  yields 

00  1 

*/ (*7?  ;M*)  -  /  e,r‘  r<~ * V*W  (17) 

-00 

where 

V4  =  k?D,{R.  ..r(l+^))  +  kiD,(R.  *o**(l-~“))  (18) 

-k^lDJp'+R.  s0^+Dp(p'-R,  soin-D^—R.  »o *)-D0'-^R,  #0r)] 

We  have  introduced  two  physical  scales  of  the  medium  which  are  defined  at  a 
wavenumber 
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These  scales  are  the  coherence  length  of  the  field,  <0, 

*oa  D„(t0)  «=1 


(20) 


and  the  size  of  the  scattering  disk,  R„ 


R. 


(21) 


The  coherence  length  is  important  in  that  the  cross-spectrum  is  negligible  for  spatial 
frequencies  k  >«o'1,  so  that  the  combination  »0k  is  always  less  than  unity.  The  scatter¬ 
ing  disk  size,  Z?,-1,  sets  the  scale  that  separates  high  and  low  spatial  frequencies.  Notice 
that  the  combination 


Rf  m  R.  »0 


(22) 


appears  in  (18).  RF  is  known  conventionally  as  the  Fresnel  radius;  it  determines  the 
intensity  spatial  scale  in  weak  scattering.  Unfortunately,  the  final  integral  in  (17)  can¬ 
not  be  done  exactly.  Our  approach  is  to  find  two  separate  approximations  for  high  and 
low  spatial  frequency.  First  we  consider  low  spatial  frequencies. 


When  jc  <l/fl„  the  0*  Fourier  transform  has  substantial  contributions  from 
everywhere  in  the  integration  plane,  and  over  most  of  this  region,  ft  ^>R,  »0k.  It  can  be 
seen  from  (18)  that  Vt  will  be  dominated  by  the  first  two  terms  on  the  right  hand  side 
of  that  equation  since  the  ff  dependent  terms  nearly  cancel.  In  this  circumstance,  the 

dominant  behaviour  is  0'  independent.  Since  the  0'  dependent  terms  in  V4  are  small, 
the  exponential  may  be  expanded  in  a  Taylor  series 


-  |  [k?  D,tf+6')+  ft?  D,(i»-P)) 


(23) 


£o7T  (Ms)’ 


where  V  =  /?,«0**and  S'  =(8k/2H)R,t0i?.  Inserting  the  above  expansion  into  (17)  yields 
explicit  expressions  for  the  terms  in  the  low-spatial-frequency  series, 

«*;(*:«  ;Ma)  -  S  */!$,(*V?;Mj)  (24) 

•  -o 
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The  structure  function  can  be  expressed  in  terms  of  a  spectrum  by 

00 

Z>„(r)  . 2 / *„(**') fl -«•**'• r]</V  (25) 

•^0 

with  as  the  power  spectrum  of  a  quantity  which  is  the  refractive  index  fluctuations 
integrated  through  the  screen.  The  two  leading  terms  in  the  low-spatial-frequency 
series  are 

(26) 

and 

=  ^exp  - 1.  [ft?  D„((1+-^R  «„*•)+*?  <o**)]  (27) 

*„(«*)  sin2(/?2/c2/2)  —  sin^-j  /?/«2/4) 

The  correction  (27)  modifies  the  behaviour  of  the  leading  order  expression,  (26),  for 
0</c<l/i?,.  Notice  that  the  Fresnel  radius  squared,  R$^R,s0  appears  in  the  sin2 
terms,  revealing  these  terms  as  the  usual  Fresnel  filter.  For  weak  scattering  (R,  <Rr ) 
the  Fresnel  filter  cuts  off  the  spectrum  at  large  spatial  frequencies.  However,  in  strong 
scattering,  the  exponential  term  in  (27)  provides  the  cutoff. 


We  now  turn  our  attention  to  high  spatial  frequencies.  It  is  convenient  to 
separate  V4  into  Vi  =  Vj+V^+Vf1,  with 

vum  ■  *,  k2[D  ,(?+*')  +  D„{P '-*')]  (28a) 

VHiU*)  =  ft?  Dtf*P}+il  D,tf-P)-2klkiDl,tf')  (28b) 

and 


V? l(*VP)  =  kxk3  -70]  (28c) 


where  7*'  <=R,t0H*  and  ?'  =(Sk/2lc)R,  a0it.  Notice  that  Vf°  is  independent  of  p'.  The 
exact  expression  for  the  cross-spectrum  is  then 


In  the  monochromatic  case  (M=0,?«0  )  at  high  spatial  frequencies  and  strong 
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scattering,  the  terms  corresponding  to  V?  dominate  the  integral  for  the  spectrum. 
[Gochelashvily  &  Shishov,  1975;  Rumsey,  1975).  Rumsey  showed  that,  over  the  impor¬ 
tant  region  of  integration,  the  ft  dependence  of  Ff*  is  negligible  compared  to  that  of 

Vj,  giving  the  intensity  spectrum  as  the  Fourier  transform  of  cxp(— *lJt3/>#<(^)).  This  is 
because  the  quantity  i0k  is  of  the  order  of  unity,  so  that  'f&R,  ,  while  the  important 
region  of  integration  is  P'&«0<R,.  When  6k  is  not  zero  (but  small),  the  high-spatial- 
frequency  behaviour  of  the  cross-spectrum  is  still  controlled  by  VJ  for  the  same  reasons. 

Note  that  these  terms  do  not  dominate  the  behaviour  of  Vt  for  all  fl'  and  1C,  but  only  for 
ft  <^R,i0k.  Unlike  the  monochromatic  case,  when  ft  is  zero  Vt  is  non-zero.  However  V4 
still  increases  with  increasing  ft,  so  that  the  Fourier  transform  integrand,  is  larg¬ 

est  when  ft  is  small.  The  approximate  domain  of  /?',  in  which  the  Fourier  integrand  is 
large,  is  where  f  <  «0.  In  this  regime  the  Vf’s  are  roughly  independent  of  f,  so  that 
the  integrand  is  dominated  by  V?.  The  limiting  value,  p&s0,  leads  to  the  requirement 

«  (30) 

wherein  we  expect  the  ft  behaviour  of  e~V*^  to  be  controlled  by  the  D^ff'+V)  +Dt,(fi'-6i) 

terms.  This  is  verified  explicitly  in  Appendix  1.  Thus  Vj1  is  small  compared  to  V? 
allowing  us  to  expand  the  final  exponential  in  a  Taylor  series  yielding  a  corresponding 
series  for  the  intensity  cross-spectrum 

*,«R ;Ma)  -  £  (31) 

•  >0 

The  leading  term  is 

00 

*,%(*, R-Kk*)  =  e~V*°  J  «<r- *  (32) 

-^0 


exp  10,(0'+^  R.  R,s0k-)}  Pft 

For  the  monochromatic  case,  the  corresponding  term  is  commonly  called  the  high- 
spatial-frequency  approximation. 

Writing  the  structure  function  in  terms  of  the  spectrum  as  in  (25),  the  second 
term  in  the  series  expansion  of  exp[--^-  Vf1]  is 

it 


.ivf> 


KUL 

2 


This  yields  the  first  correction  to  the  leading  approximation,  (32), 


♦aw**)  **•  •""> 

*  -00 


exp  -^[D'(fi'+^R.s0l?)+Dp(?'-yj:R..0*))  »n2  d'fi'J**' 

where  V*0  is  given  in  (28b).  This  correction  term  is  important  in  determining  the  con¬ 
vergence  properties  of  the  series  (31)  for  high  spatial  frequencies. 


4.  Comparison  to  the  Gaussian-Field  Approximation 


In  strong  scattering,  the  coherence  length  of  the  field,  «0)  is  much  smaller  than 
the  size  of  the  scattering  disk.  This  implies  that  the  field  at  the  observation  point  is 
the  sum  of  very  many  contributions  from  the  scattering  disk.  If  these  contributions 
were  independent,  then  by  the  central  limit  theorem  the  field  would  obey  Gaussian 
statistics.  If  the  field  is  a  zero-mean  Gaussian  random  process  then  the  intensity  corre¬ 
lation  function  is 


Car(*i& s>fl;Ma)  =  <  I{x\,R  ;*i)>  <  /(x^fi;*2)>  +  |r2(x*1,xy?;*1,*2)|2  (35) 

In  the  thin  screen  model  and  an  incident  plane  wave,  the  mean  intensities  in  the  first 
term  are  independent  of  position  and,  by  definition,  are  unity.  Because  of  translation 

invariance  the  second  term  is  only  a  function  of  0  =  x*|-S£.  Inserting  r2  from  (12)  and 
transforming  over  0  yields  the  GF  approximation  to  the  intensity  cross-spectrum 


=  *(*')  + 


)2<|>2 


exp  - 


We  wish  to  compare  this  GF  result  with  our  previously  derived  approximations: 
(26)  4-  (27)  for  small  k,  (32)  +  (34)  for  large  k.  We  see  that  (36)  resembles  (26)  +  (32) 
except  that  the  factor  exp[-Vf°]  has  been  replaced  by  the  phase-shift  factor 
exp(-(M /2F)9*8].  However  there  is  a  problem  with  the  low-spatial-frequency  behaviour 
of  the  GF  approximation.  On  quite  general  grounds  it  can  b’e  argued  that,  for  low  spa¬ 
tial  frequencies,  the  cross-spectrum  consists  of  a  delta  function  piece  plus  a  term  that 
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must  vanish  as  the  spatial  wavenumber  goes  to  zero  [Tatarskii,  1971).  If  (36)  is  applied 
for  small  x,  the  delta  function  is  correct  but  the  second  term  goes  to  a  non-zero  con¬ 
stant  for  small  x.  Compare  this  situation  with  our  results  (26)  and  (27).  The  result  (26) 
correctly  gives  the  delta  function  at  the  origin.  Our  correction  (27)  substantially 
modifies  the  spectrum  for  non-zero  x,  but  satisfies  the  requirement  of  vanishing  as 
k  — *  0.  For  a  power-law  medium,  such  that  4>M«x-®,  with  or  <4,  it  is  seen  that  (27)  van¬ 
ishes  as  x4_<*. 

To  evaluate  the  GF  approximation  at  high  spatial  frequencies,  we  are  interested 
in  whether  exp|~($Jfc/2F)2<l>2}  is  a  good  approximation  to  exp[-Vf0].  Consider  a  pure 
power-law  phase  structure  function 


•  <  ^ otter 


(37a) 


*02z>„(r)=  <t>2  «>/„,„ 


(37b) 


where  «0  *s  the  coherence  length  of  the  field,  and  ltnUr  is  the  outer  scale  of  the  spectrum 
of  phase  irregularities.  In  most  cases  of  practical  interest,  the  outer  scale  is  much 
larger  than  the  size  of  the  scattering  disk. 

Since  the  phase  structure  function  saturates  to  4>2  for  scales  larger  than  the 
outer  scale,  (32)  gives  the  same  result  as  (36)  only  for  spatial  frequencies  such  that 
x  >(/,„«r/7?.)«6'1-  The  only  non-negligible  portion  of  the  spectrum  is  for  x  <!/*„,  leading 
to  the  conclusion  that  the  GF  approximation  is  only  valid  for  a  significant  portion  of 
the  spectrum  when  the  outer  scale,  f„«r,  is  small  compared  to  the  size  of  the  scattering 
disk.  However,  in  virtually  every  case  of  practical  interest,  the  scattering  disk  is  small 
relative  to  any  estimates  of  an  outer  scale.  In  these  cases,  the  GF  approximation 
becomes  valid  only  after  the  cross-spectrum  has  dropped  to  a  negligible  value. 


5.  Comparison  to  Neglecting  the  Phase-Shift  Factor 


As  mentioned  previously,  various  arguments  have  been  given  for  using  (36)  for 
the  cross-spectrum  but  neglecting  the  phase  shift  factor, 

OO  i  * 

♦or(«:*;Ma)  «  /  R>  R>  *o**)l  W  (38) 

Dashen  (1979)  considered  this  problem  theoretically  for  both  a  phase  screen  and  for  an 
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extendcd  random  medium  and  presented  results  valid  for  infinitesimal  frequency 
differences.  There  are  two  different  regimes  in  strong  scattering:  "  partial”  and  M  full" 
saturation.  When  the  size  of  the  scattering  disk  is  larger  than  the  outer  scale  of 
medium  fluctuations  (full  saturation)  the  contributions  are  independent  and  the  GF 
approximation  is  valid  to  first  order.  However,  when  the  scattering  disk  is  smaller  than 
the  outer  scale  (partial  saturation),  fluctuations  with  a  size  larger  than  that  of  the  disk 
contribute  a  random,  coherent  phase  to  the  field  so  the  many  contributions  are  not 
independent.  Since  this  coherency  only  affects  the  rms  phase,  this  lack  of  independence 
will  not  affect  any  monochromatic  intensity  statistic  but  will  be  important  for  multifre¬ 
quency  statistics.  For  a  medium  characterized  by  a  fluctuation  spectrum  with  a  spec¬ 
tral  index  less  than  four,  the  small-scale  fluctuations  can  cause  large  (saturated)  inten¬ 
sity  fluctuations,  but  the  rms  phase  is  dominated  by  the  large-scale  fluctuations.  For 
small  enough  frequency  differences,  this  phase  cannot  affect  the  intensity,  providing  a 
first-order  justification  for  dropping  the  exponential  factor  containing  the  rms  phase. 
Further  experimental  confirmation  of  this  result  was  provided  from  ocean-acoustic  data 
[Flatte',  1983]. 

Our  more  accurate  expression  (32)  implies  that  neglecting  the  phase-shift  factor 
is  valid  over  the  portion  of  the  spectrum  for  which  Vj0  Cl.  Using  the  largest 
significant  k  value  (a^1 )  in  (28b)  leads  to  the  requirement 

*?I)M(H.(l+^))+*?f)|,(/?.(l-||))-2*I*2^(/?.)  <1  (39) 

For  a  power  law  structure  function  this  expression  can  be  simplified  by  performing  a 
small  6k  expansion  (  valid  for  6k  <2F  )  leading  to 

(^U  Cl  (40) 

where  U  is  defined  by 

U  S  k$D,(R.)  &  '  (41) 

•° 

This  requirement  implies  that,  as  the  scattering  becomes  stronger,  the  region  (in  6k  )  for 
which  (38)  is  valid  becomes  smaller.  This  is  not  a  severe  restriction  because  the  fre¬ 
quency  decorrelation  bandwidth  decreases  as  the  scattering  becomes  stronger,  limiting 
the  values  of  6k  which  enter  (40).  We  set  this  limiting  value  to  be  the  decorrelation 
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bandwidth,  which  we  estimate  by  looking  at  the  asymptotic  (large  k  )  form  of  (38). 

For  large  k,  there  are  two  important  regions  of  integration  in  (38):  near 

where  the  integrand  is  largest,  and  near  =  ±(6k/2k)R,»^t,  where  one  of  the  structure 
functions  vanish.  In  Appendix  2,  it  is  shown  that  the  second  region  dominates  the 
asymptotic  behaviour,  leading  to 


exp— 


Ms 

2 


D„( 


oo  ^1^2  _i> 

^-R.s o/r)cos  J  (42) 

f  2Jc  I  -OO 


The  integral  term  in  (41)  is  just  the  angular  spectrum  and  its  important  scale  is  e0.  The 
cosine  term  oscillates  rapidly  for  k  =  V 6k/2Tc  / RF .  The  frequency  decorrelation 

bandwidth  can  be  estimated  as  that  6k  for  which  Using  (41)  and  (21)  this  leads 

to 


6kh  «  jrFl/->/r  (43) 

Inserting  this  bandwidth  in  (40)  yields  the  criterion  for  neglecting  the  phase-shift  factor 

t/i-8/r<l  (44) 

Thus,  for  p  <2,  the  approximation  of  neglecting  the  mean-square  phase  shift  improves 
as  the  strength  of  scattering  increases. 

It  has  been  shown  that  the  contribution  to  the  scintillation  index  from  the  high- 
spatial-frequency  correction  term  is  proportional  to  Ux~*l*  [Prokhorov  et  al,  1975). 
Thus  if  the  correction  term  to  the  scintillation  index  is  small,  it  is  appropriate  to  ignore 
the  phase-shift  factor  over  the  main  portion  of  the  spectrum.  On  the  other  hand,  if 
correction  terms,  such  as  (34),  are  not  neglected,  then  it  is  inappropriate  to  neglect  the 
phase-shift  factor,  since  it  generates  corrections  of  the  same  order-of-magnitude. 


6.  Summary 

In  the  previous  sections  we  derived  two  different  series  for  the  intensity  cross¬ 
spectrum:  one  describes  the  low-frequency  behaviour  while  the  other  describes  the  high- 
frequency  behaviour.  Keeping  the  first  two  terms  in  each  series  gives  expressions  which 
approximate  the  cross-spectrum  in  the  two  respective  regions.  The  relevant  formula  are 
given  in  (26), (27), (32),  and  (34).  The  phase-shift  factor  in  (32)  containing  V$°  correctly 
describes  the  effects  of  any  coherent  contributions  from  the  scattering  disk  for  all  values 


of  the  outer  scale  and  for  arbitrary  frequency  differences. 

The  GF  approximation  should  not  be  used  when  the  outer  scale  is  much  larger 
than  the  diameter  of  the  scattering  disk,  and  the  reduced  approximation,  (40),  should 
only  be  used  when  1.  In  very  strong  scattering  (large  U),  when  the  scintillation 

index  has  nearly  saturated  to  unity,  (40)  provides  a  good  approximation  over  the  main 
portion  of  the  spectrum. 


Acknowledgements:  This  research  was  funded  by  DARPA. 


References 


Azar,  Z.,  H.M.  Loebenstein,  G.  Appelbaum,  E.  Azoulay,  U.  Halavee,  M.  Tamir,  and  M. 
Tur,  Aperture  averaging  of  the  two-wavelength  intensity  covariance  funtion  in  atmos¬ 
pheric  turbulence,  Applied  Optics,  24,  2401-2407,  1985. 

Codona,  J.L.,  D.B.  Creamer,  S.M.  Flatte'  ,  R.G.  Frehlich,  and  F.  Henyey,  Solution  for 
the  fourth  moment  of  waves  propagating  in  random  media,  submitted  for  publication, 
1985. 


Cole,  T.W.  and  O.B.  Slee,  Spectra  of  interplanetary  scintillation,  Nature,  285,  93-95, 
(1980). 

Cordes,  J.M.,  and  J.M.  Weisberg,  and  V.  BoriakofT,  Small-scale  electron  density  tur¬ 
bulence  in  the  interstellar  medium,  Ap.  /,  288,  221-227,  1985. 

Dashen,  Roger,  Path  integrals  for  waves  in  random  media,  J.  Math.  Phys.  20,  894-920, 
1979. 

Flatte'  ,  S.M.,  Wave  propagation  through  random  media:  Contributions  from  ocean 
acoustics,  Proc.  IEEE,  71,  1267-1294,  1983. 

Ginzburg,  V.L.,  and  L.M.  Erukhimov,  On  the  influence  of  electronic  concentration 
fluctuations  of  the  dispersion  measure  and  rotation  measure,  Astrophys.  Space  Sci.,  11, 
362-372,  1971. 

Gochelashvily,  K.S.,  and  V.I.  Shishov,  Saturation  of  laser  irradiance  fluctuations  beyond 
a  turbulent  layer,  Optical  and  Quantum  Electronics,  7,  524-536,  1975. 

Gurvich,  A.S.,  and  V.  Kan,  Fluctuation  of  intensity  at  two  wavelengths  in  a  turbulent 
medium,  Izv.  Vysh.  Ucheb.  Zaved.  Radiofiz.,  22,  843-847,  1979.  English,  Radiophysics 
and  Quantum  Electronics,  22,  585-587,  1979. 


Gurvich,  A.S.,  V.  Kan,  and  VI.  V.  Pokasov,  Two-frequency  fluctuations  of  light  inten¬ 
sity  in  a  turbulent  medium,  Optica  Acta,  2d,  555-562,  1979. 

Jakeman,  E.,  G.  Parry,  E.R.  Pike,  and  P.N.  Pusey,  The  twinkling  of  stars,  Contemp. 
Phys.,  19,  127-145,  1978. 

Lee,  L.  C.,  Strong  scintillations  in  astrophysics.  IV.  Cross-correlation  between  different 
frequencies  and  finite  bandwidth  effects,  Ap.  J.  206,  744-752,  1976. 

Popov,  M.V.  and  V.A.  Soglasnov,  Observations  of  interstellar  scattering  of  radio  emis¬ 
sion  of  pulsar  PSR  0329454  at  102.5  MHz,  Astron.  Zh.,  61,  727-732,  1984.  English,  Sov. 
Aslron.,  28,  424-427,  1984. 

Prokhorov,  A.M.,  F.V.  Bunkin,  K.S.  Gochelashvily,  and  V.I.  Shishov,  Laser  irradiance 
propagation  in  turbulent  media,  Proc.  IEEE,  63,  790-811,  1975. 

Rickett,  B.J.  Frequency  structure  of  pulsar  intensity  variations,  Nature,  221,  158-159, 
1969. 

Rumsey,  V.H.,  Scintillations  due  to  a  concentrated  layer  with  a  power-law  turbulence 
spectrum,  Radio  Science,  14,  107-114,  1975. 

Shishov,  V.  I.,  Frequency  correlation  of  scintillations,  Izv.  Vysh.  Ucheb.  Zaved.  Radio fiz. 
18,  423-433  (1973).  English,  Radiophysics  and  Quantum  Electronics  16,  319-327,  1973. 

Tatarski,  V.  I.,  The  Effects  of  the  Turbulent  Atmosphere  on  Wave  Propagation, 
National  Technical  Information  Service,  Springfield,  Va.,  1971. 


9— 18— '85 


Appendix  1 


We  show  that  f?  dependence  of  Vf 1  is  negligible  compared  to  that  of  Vj  where, 
from  (28), 

V? (Ala) 

and 

Vf*  =M2[2  Djff-Djfl+fi-Djfi-?)]  (Alb) 


where  i=R,s q**  and  6=(6k/2k)R,srf?.  We  have  argued  in  section  3  that  the  dominant 
region  of  integration  is  for  small  0'  •</?,  «0/c.  In  that  case,  the  two  terms  in  (Alb)  may 
be  expanded  as 

^(®*)lr-«,»0r  +  (A2) 

The  second  term  may  be  bounded  by 

<  0l2\DM"(R,  «0«)|  (A3) 


?jz 


DJF)  I 


r-R,.0r 


Dh{R,  »o*'+Pl)+Dlt(R,  a0K*-p')  =  2 Dm(R,  «0k*)+ 


giving  the  leading  01  dependence  of  V* 1  as 

V?  »  M2/?2  D„"{R,  »0k)  (A4) 

In  verifying  that  Vj  controls  the  fP  dependence  we  will  examine  three  cases 


Case  1 : 

^zR,»ok  <0* 

2k 

Case  2: 

01  <  —~R,  sqk 

2k 

(A5) 

Case  3: 

j^RssoK  «  p 

2k 

In  the  first  case,  0  is  large  compared  to  SkR,  e0n/2k  giving  the  leading  dependence 
of  Vj  as 

V?  «  2klk3D„{0>)  (A6) 

so  that  we  require 

(A7) 
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Using  a.  model  power-law  structure  function,  Jfc  *£>,,(«)  =(*/<0)',  this  requirement  becomes 


< 


2 

P(P-1) 


(AS) 


Since  p  is  typically  in  the  range  0  <p  <  2,  we  find  an  upper  bound  for  the  left-hand  side 
by  setting  p1  =  «0  and  k  =«^_1  resulting  in  the  requirement 


-E-ii — ^  yi2-p)/p 
2 


(A9) 


where  U  is  the  strength  of  scattering  parameter,  U  =  k$  D,,{R,),  with  Jfc0  defined  in  (19). 

Therefore,  the  term  controls  the  p  dependence  when  the  strength  of  scattering 
parameter,  U,  is  large.  When  p  is  1  or  2,  our  argument  breaks  down.  For  p  =1  the 
approximation  still  holds  by  another  argument,  while  for  p  =2,  there  are  no  scintilla¬ 
tions. 


In  the  second  case,  P'  is  small  compared  to  6kR,  *0k/2F,  so  that  a  Taylor  expansion 
of  the  V°  term  yields  the  leading  p'  dependence 

V}  «  k^D/i^R.  *0k)  (A10) 

This  leads  to  the  requirement 

D;\R. «,*)  CDSl^R.  80k)  (All) 

which  for  a  power-law  structure  function  is  equivalent  to 

P-2 

(A12) 

Since  p  is  less  than  2,  this  implies  that  6k  is  small  compared  to  twice  the  mean  of  the 
wavenumbers. 

Finally,  in  the  third  case,  p1  is  about  the  same  size  as  6kR,  «0k/2F  so  that  the 
requirement  for  neglecting  the  p'  dependence  of  V$  becomes 

» 

D„”  <&D^R,8 0(c)  (A  13) 

which,  for  a  power-law  structure  function,  leads  to 


6k  D 


1  « 


6k_ 

23F 


(w/f)  <7(^  (AM) 

and  is  easily  satisfied  for  small  wavenumber  differences. 


In  all  three  case  considered  above,  the  combination  of  structure  functions  in  VJ 
dominate  the  ft  dependence  if  U  is  large  and  6k/2t  is  small. 


i 

i 

i 

I 

i 
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Appendix  2 


As  mentioned  in  section  4,  in  (38)  there  are  two  important  regions  of  integration 
for  large  k.  One  is  for  small  0  where  the  integrand  is  largest.  The  other  is  where  one  of 

the  structure  functions  vanish,  0=±(Sk/2T)Rt»OK.  Near  0=Q  the  exponent  in  the 
integrand,  Vf  may  be  approximated 


2z>„(-^?,ao/n+ 


da 


(Bl) 


which  leads  to  rapidly  decaying  component  of  the  cross-spectrum.  For  the  other  region 

of  integration  we  approximate  the  integrand  of  (37)  in  the  vicinity  of  ff' =(6k /2X)R,»0k* 
to  be  dominated  by  the  “fast”  variation  of  the  structure  function  near  its  zero  and  the 

“slow”  variation  of  the  other  portion  of  the  exponent.  Including  the  /T=0  contribution 
gives  the  cross-spectrum  asymptotically  as 


***»  D  r J*  R  &  <» 
~D+i r  V1 


exp 


!^(0j  rH*/2)F)«,.0r] 


(B2) 


J^Zd&k,  m) 

+2  e  *  s  0  cos 


Hr}** 

2JF  F 


-oo 

For  a  power-law  structure  function,  the  first  term  falls  off  like 


CxK^t~C^~7  (B3) 

where  C,  and  C%  are  constants.  The  integral  in  the  second  term  is  just  the  angular 
spectrum,  which  for  p  <2,  falls  off  like  k -v-*.  Since  this  is  a  much  slower  fall  off  than 
the  0-0  contribution,  the  asymptotic  behaviour  of  the  cross-spectrum  is  dominated  by 

the  points,  0=±{8kf2t)Rfi?. 
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Abstract 

Differential  equations  for  all  moments  of  the  field  of  a  wave  propagating  through  a 
random  medium  are  derived  under  the  parabolic  approximation  and  the  Markov 
approximation,  hut  including  anisotropy  in  the  random  medium  and  a  deterministic 
background  refractive  index.  Mathematical  equivalence  is  demonstrated  between  these 
moment  equations  and  path-integral  expressions  for  the  moments  obtained  under  the 
same  approximations.  A  discussion  of  approximations  that  are  weaker  than  Markov  is 


given. 


L  Introduction 

Many  problems  in  wave  propagation  through  random  media  concern  phenomena  in 
which  there  is  no  significant  backscatter,  so  that  a  parabolic  approximation  may  be 
made  to  the  wave  equation. ^  In  these  cases  a  further  approximation,  called  the  Markov 
approximation,!8!  leads  to  relatively  tractable  mathematical  expressions  for  moments  of 
the  field  that  can  be  used  for  practical  calculations.  Two  quite  different  formalisms 
have  been  used  in  this  context:  the  moment-equation  and  path-integral  techniques. 

A  path-integral  expression  for  a  general  moment  of  the  field  of  a  wave  propagating 
through  an  inhomogeneous,  anisotropic  medium  in  the  presence  of  a  deterministic  back¬ 
ground  refractive  index  has  been  derived,!3!  and  the  expression  has  been  used  for  specific 
calculations.!4,5,8! 

Moment  equations  in  coordinate  representation  have  been  derived  for  homogene¬ 
ous  isotropic  media  in  the  absence  of  a  deterministic  background.!3!  Treatments  of  inho¬ 
mogeneity,  anisotropy,  and  deterministic  background  by  moment-equation  techniques 
have  heretofore  been  confined  to  special  cases  involving  the  first  and  second 
moments.!7,8! 

We  present  here  general  moment  equations  in  coordinate  representation  that 
account  for  inhomogeneity,  anisotropy,  and  deterministic  background,  but  require  the 
Markov  approximation.  We  derive  these  equations  using  the  time-ordered-product 
method  of  Van  Kampen,!8!  which  also  provides  a  derivation  of  equations  that  are  valid 
under  conditions  more  general  than  the  Markov  approximation.  The  modified  equations 
are  more  complicated  than  those  that  require  the  Markov  approximation:  a  special  case 
was  previously  derived  by  Besieris  and  Tappert.!10! 

We  also  show  that  our  new  genera]  moment  equations  derived  under  the  Markov 
approximation  are  mathematically  equivalent  to  the  path-integral  expressions  for  the 
moments  that  have  been  previously  presented.  Thus,  the  two  popular  formalisms, 
under  the  Markov  approximation,  are  not  different  in  content. 
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The  plan  of  the  paper  is  as  follows:  in  Section  U  we  establish  notation,  present  our 
new  moment  equations,  and  present  path-integral  expressions  for  the  moments  in  simi¬ 
lar  notation.  In  Section  III  we  establish  the  mathematical  equivalence  between  the  two 
techniques.  In  Section  IV  we  present  the  derivation  of  our  moment  equations,  and, 
along  the  way,  derive  the  modified  equations.  In  Section  V,  for  completeness,  we 
rederive  the  path-integral  expressions  for  the  moments.  In  Section  VI  we  comment  on 
the  use  of  different  coordinate  systems  (such  as  cylindrical  or  spherical)  in  the  writing  of 
moment  equations.  A  summary  concludes  the  paper. 


11.  Notation  and  Markov-Approxlmatlon  Results 

Consider  waves  (ravelling  predominantly  in  the  »  direction.  Let  f*  be  a  transverse 
coordinate  (e.g.  two-dimensional,  but  in  fact  general),  and  *  be  a  reference  wave 
number  (Jr  ■»  2tu/C0),  where  u>  is  the  wave  frequency  and  C0  is  a  reference  wave  speed). 
Express  the  full  wave  field  as 

«(F,*,I)  -  i&(r,i)exp[tfc(* -<V)]  (1) 

Let  the  wave  speed  (a  function  of  position  only)  be 

C(?,x)  =  C0  [  1  - 2 1/0(**)  ~ *) ]  14  »  C0[l  +  t/0(r)+/i(r,e)]  (2) 

where  £/„  represents  the  deterministic  background  and  p  represents  the  fluctuating  ran¬ 
dom  medium,  assumed  to  be  a  realization  of  a  zero-mean  Gaussian  process. 

Then,  the  parabolic  equation  (in  rectangular  coordinates)  for  the  reduced  wave 
function  is: 

ikd,i>  .  -  ±  W  +  *5C0(*>  +  z)i>  (3) 

where  V2  is  the  transverse  Laplacian. 

A  moment  T  is  the  ensemble  expectation  value  of  a  product  of  i&’s  and  i>  *’s  where 

each  ^  or  is  evaluated  at  a  different  position  fj  and  wavenumber  kf.  We  write,  in 
abbreviated  form, 

r„.  »  <*?■■  H>m  +n+ 1  •  *  *  >  (4) 

Define  an  operator  L0  such  that 

i.  -  (  -jVj+t/iJ  (S) 

The  terms  that  apply  to  the  i>'s  use  the  plus  sign  and  those  that  apply  to  the  i>*'a  use 
the  minus  sign.  The  subscript  /  requires  that  V*  operate  only  on  S}  and  m  U0{§}). 

Define  the  important  combination  of  fluctuation  quantities  as 

M(*)  -  ±  kjp(s;.,z)  (6) 
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Our  general  moment  equation  under  the  Markov  approximation  can  be  written 


».r„M  -  ]  u  <M(.)M«(.')>r„(.)  (7) 

—  00 

where  is  obtained  by  evaluating  M(t)  with  all  the  IJ  at  t  shifted  by  the 

transverse  distance  that  a  deterministic  ray  through  moves  in  travelling  from  x 
to  x *  (see  Figure  1).  In  other  words  is  evaluated  at  point  B:  i.e.  gj  -  rrtt  {x1) 

where  the  ray  is  forced  to  go  through  f?(«).  The  particular  ray  is  determined  not  only 
by  the  local  position  (ij  ,z),  but  also  by  the  initial  conditions  on  the  moment;  for  exam¬ 
ple,  the  location  of  a  point  source,  or  the  direction  of  a  plane  waveJ11)  The  unphysical 
assumption  of  delta-correlated  medium  fluctuations  along  the  propagation  direction 
would  imply  that  would  be  evaluated  at  point  C:  i.e.  z?(z)  (and  x1).  In  the  iso¬ 

tropic  case  (or  in  the  case  of  propagation  along  a  principal  axis  of  the  anisotropy)  the 
difference  between  evaluating  M .*#(**)  at  ay(z)  and  (z*)  is  negligible,  and  the  delta- 
correlated  assumption  is  adequate.  In  the  anisotropic  case,  the  necessity  of  defining  the 
unperturbed  ray  makes  (7)  somewhat  complicated  to  apply  for  general  initial  condi¬ 
tions.  However,  since  (7)  is  a  linear  equation,  superposition  can  be  used  whether  the 
source  is  a  point,  an  incident  plane  wave,  or  an  arbitrary  coherent  or  incoherent  sum  of 
point  sources.  Equation  (7),  which  is  one  of  the  principal  results  of  this  paper,  is 
derived  in  Section  IV. 

We  now  turn  to  the  path  integral  method.  Equation  (3)  has  the  formal  solution 

0=/£»r(z)  c*  (8) 

where  /  Di’(x)  means  integration  over  paths,  f*(z)  is  a  transverse  vector  indicating  the 
position  of  the  path  at  z,  and 

*  -  */*  |  (tt) (») 

In  order  to  obtain  a  given  moment,  expressions  like  (8)  (or  its  complex  conjugate)  are 
multiplied  together,  and  the  ensemble  average  is  taken: 

m±m  H*iS] 

r-  -jna*W<«'  > 

/-I 


(10) 


The  Markov  approximation  yields  (See  Section  V): 


M> 

•  —00 

(11) 

We  show  in  the  next  section  that  the  moment  equations  (7)  and  the  path  integral 
expressions  (l  1)  are  mathematically  equivalent. 
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Equivalence  of  Path  Integral  and  Moment  Equations  under  the  Markov 
Approximation 


We  follow  the  technique  that  Feynman!12*  used  to  show  that  his  path-integral 
expression  for  nonrelativistic  quantum  mechanics  is  equivalent  to  the  Schrodinger  equa¬ 
tion.  The  key  to  this  demonstration  is  an  understanding  of  how  the  important  paths 
behave  transversely  as  they  move  in  *  from  a  particular  point.  Feynman  found  that 
these  paths  resembled  random  walks  in  that 

i*-v)-r(*)i  ~  (*'-*)*  (12) 

as  x '  gets  close  to  x.  Given  this  behavior,  it  is  easy  to  expand  (11)  in  a  Taylor  series 
and  obtain  a  differential  equation  which  will  turn  out  to  be  (7).  We  give  the  demonstra¬ 
tion  of  (12)  in  the  Appendix. 

The  path  integral  is  defined  as  the  limit  of  an  integration  over  a  set  of  “phase 
screens.”  These  screens  are  at  values  xN  =  Nix.  The  derivative  at  x  =  xN  is 

&t 

defined  as  (f*(zw  +  ix)  -  z*(xN)) /Sx  =  The  limit  is  taken  after  the  integrals 

are  evaluated.  The  differential  equation  is  obtained  by  considering  the  integral  over  the 
very  last  phase  screen.  The  last  integral  in  (11)  can  be  written  in  terms  of 
xy  {R  —  Sx)  and  Sj  »=*  f/(J?).  Also,  we  define  Bz}  s  x}  —  £}'.  Then  r„,  can  be 
expressed  as: 


r  ..({DM  = 


N/n^'exP 
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-jfdx'<M(R)M,w(*')> 


rmMn,R-ix) 


(13) 


where  {f*}  denotes  the  set  of  m+n  Sj’s.  The  first  term  in  the  exponent,  ±~^~  -1 —  , 

2  ox 

is  0(1)  for  small  6x,  because  of  (12).  The  exponent  of  the  remaining  terms  can  be 
expanded,  since  they  have  an  explicit  6x,  as  well  as  higher  order  terms.  This  results  in 


i 
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r«.«r>,*) 


(H) 


n/  n  ('«*}  exp  {±  ik,m  )V**}) 


l-Sz  S±«’*/l/oy(fi)  +  “-  /  iz'<M(R)Mtkifl(zf)>  r„.«r'}  ,* -*2) +0(5,») 
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We  now  have  a  relationship  between  the  moment  at  R  and  the  moment  at  R  —  6z, 
which  we  derived  from  our  path-integral  expression.  But  since  the  moment  is  a 
differentiable  function  we  can  find  another  relationship  by  Taytor  expansion  as  follows: 

=  i-szdR  -  ’Zszj'Vj  +  y(5>*rv/)8]  r«..({  **).*)  (15) 


+  0(5i 2  ) 


Substituting  (15)  into  (14)  we  find 

I\„({  **},*)  «N/n(rf5i}  exp{±.*J(5*JS)/252}) 


-52  E  ±*iUoil*)  +  j  !j*'  <  M(R)MMfl(z')  >  |r(„({r}>/?)  +o(5,2) 


The  term  linear  in  E is  odd  in  5 ij  and  therefore  gives  zero  due  to  the  5fj 

integral.  The  term  that  is  quadratic  in  5*^  can  be  integrated  by  parts,  yielding,  to 
order  52: 


r-«r>,s)  -  N/n  ('*»;«xp{±iMJr,)72J«}) 


-{  J  /.'<*/(*)«* „M> 
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r„. an.R)  . 


(17) 


The  only  way  (17)  can  be  true  for  all  Sz  is  for  tbe  coefficient  of  6z  within  the  curly 
brackets  operating  on  T,,,  to  give  zero.  Therefore,  setting  R— *z, 


».rM((r},«)  -  +  rm.({r},*) 


7  /  ^,<M(i)M,ikiyi(zO>rB.({r},o 

*  -oo 


(18) 


which  is  identical  to  (7),  as  required.  Thus,  we  have  derived  the  moment  equation  (7) 
from  the  path  integral  expression  (11).  This  shows  that  the  path-integral  expression 
(11)  is  a  solution  of  the  moment  equation  (7)  and  bonce  the  two  techniques  are 
equivalent. 


IV.  Moment-Equation  Derivation 

We  derive  our  moment  equations  by  the  method  of  Van  Kampen.M  The  advantage 
of  his  method  is  that  the  physical  basis  for  each  approximation  is  readily  apparent.  He 
bases  his  method  on  techniques  that  were  developed  for  quantum  mechanics. 

We  shall  find  that  the  Markov  approximation  requires  that  the  dimensionless 
number  be  small  where  L,  is  the  medium  correlation  length  in  the  direction  of  the 
wave  propagation,  and  Mt  is  the  “typical”  value  of  M,  defined  by  (6)  and  called  the 
“interaction  strength.”  For  the  first  moment  M  =  kfi,  but  for  higher  moments  M  is  the 
sum  and  difference  of  a  number  of  Jr/»* s  at  different  positions,  and  with  different  values 
of  k. 

We  start  with  the  parabolic  wave  equation  (3)  and  the  definition  of  L0  and  M,  and 
write: 

*  d.'t’*  ta*  ■  i>m  +.  =  (L0  +  i>*  ■  •  •  V-m  -w  (19) 

The  “interaction  representation”  is  defined  by: 

(tf  iW'  * '  *m  +,)/  «  *  L°‘  t*  ■  ■  ■  I’m  (20) 

and 

M,(z)  -  (21) 

With  these  definitions,  (19)  becomes 

•  i>a  ‘ ‘  1>m  +.)/  *  Af/(* W*  +*•••  +»  +.)/  (22) 

This  equation  is  linear  and  has  the  formal  solution 

{1>*  ‘  1>m  +«)/  =  T exp  (—if  Mf(z')  dz ')  T„(0)  (23) 

o 

r„,(0)  b  the  initial  condition.  The  “time-ordering”  symbol  T  requires  explanation.  One 
notices  that  Mt  b  an  operator,  not  just  a  function  of  space.  Mt{zx)  and  Af;(z3)  do  not,  in 
general,  commute.  If  they  did  the  solution  of  (22)  would  be  given  by  (23)  without  the  T 

symbol.  The  T  symbol  means  that  a  product  of  operators  to  the  right  b  not  applied  in 

the  usual  order,  but  in  such  a  way  that  operators  with  smaller  values  of  z 1  are  to  be 
applied  first.  Thus  there  b  an  ordering  in  *.  (The  T-symboi  was  invented  for  solving 
problems  in  quantum  mechanics  where  the  analog  of  the  longitudinal  direction  z  b  tbe 
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time.)  For  example, 

Ttx p(  -  •  /  is1)  - 


[Texp(  -  i  f  U,(m 0  *')]  [Texp{  -  .•  /  A/,(*0  **)] 


for  0  <  x,  *.  Another  example  is 

/W'f  (25) 

where  the  integration  region  on  the  right  side  of  (25)  which  is  k\ 

times  smaller  than  that  of  the  left  side,  cancelling  the  factor  of  k\.  Using  either  (24)  or 
(25),  one  readily  checks  that  (23)  is  a  formal  solution  of  (22). 

We  are  assuming  that  M  is  a  Gaussian  process.  The  result  that  the  expectation  of 
the  exponential  of  a  zero-mean  Gaussian  random  variable  is  the  exponential  of  half  the 
variance  follows  from  combinatorial  factors  and  remains  true  for  a  time-ordered 
exponential.  Thus 


(r„),  -  Ttx p 


-J<  /  j 


>  r„<o) 


Although  this  is  a  formal  expression  for  rm,  it  is  not  immediately  useful  for  calcula¬ 
tions,  since  there  is  no  simple  algorithm  for  evaluating  a  time-ordered  exponential  (in 
contrast  to  a  normal  exponential).  Van  Kampen  proceeds  by  differentiating  (26): 


«,( r„.),  - 


-T<M,(t)]  dxlM,{>l)>  exp  -|<  /*#,(**) 


yg 

)m ,(*»)*•  >  r„.( 
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The  M/(t)  has  the  largest  z,  so  it  is  written  in  the  proper  ordered  position.  The 
that  it  is  correlated  with,  however,  might  occur  anywhere  relative  to  the  M,(i  ")’s  in  the 
exponential.  If  LfM?  <  1,  very  little  error  is  made  by  assuming  that  the  first  two  A/,’s 
are  in  the  proper  order,  so  that  the  T  symbol  can  be  brought  through  the  first  expecta¬ 
tion  value,  yielding: 


».(r~);  -  -<«,(.)/*' M/M  xr„),  (») 

This  may  be  shown  by  expanding  the  exponential  operators  in  (26)  or  (27)  and  dis> 
cussing  the  order  of  M* s  in  each  term.  The  Nth  term  in  the  expansion  has  2 N 
occurrences  of  Mlt  and  is  of  a  magnitude 

<(/ M,  it,f>N/N\  (29) 

o 

where  typical  eigenvalues  of  the  operators  are  implied.  The  terms  beyond 

N  *  4<(/M,  *»)*  >  (30) 

o 

become  negligible  compared  to  the  original  exponential  in  (26),  so  we  have  to  deal  with 
at  most  N  pairs  of  Af/’s  from  source  to  range  i.  The  two  M;’s  in  a  correlated  pair  must 
be  within  Lf  of  each  other  to  give  a  nonzero  correlation.  The  number  of  pairs  may  be 
estimated  as 

N  »  4L9zM?  (31) 

where  Af?  is  a  typical  value  of  Mf.  (See  Figure  2  for  a  schematic  representation.)  Our 
approximation  reduces  to  saying  it  is  unlikely  to  find  a  third  occurrence  of  an  Mt  in 
between  a  pair  that  are  within  L,  of  each  other.  This  probability  is  roughly 

Probability  «  — ^ —  «  L’AfJ*  (32) 

Thus  if  the  fluctuations  are  weak  enough  (LpMt  <L  l),  the  approximation  is  valid,  and 
(28)  is  justified. 

We  call  (28)  "first  order  perturbation  theory.”  In  typical  situations,  *  is  much 
larger  than  L/>,  and  the  lower  limit  can  be  replaced  by  -oo,  making  the  equation 
independent  of  the  source  position.  Moreover,  the  integral  from  —  oo  to  a  can  be 
replaced  by  half  the  integral  from  -oo  to  oo,  when  the  correlation  is  a  much  slower 

function  of  -^(*  +  r')  than  of  *  -  a*.  The  result  is  used,  not  in  the  interaction  represent 

ration,  but  in  the  original  representation.  The  exponentials  of  (20)  and  (21)  are 
removed,  giving 


-a,r„.(«)-  /  u<  <«(.).  -'**  -*»  >  r„,(i) 

-00 

For  the  second  moment,  this  equation  is  related  to  an  expression  of  Besieris  and  Tap- 
pertJ10*  Although  their  work  was  for  the  second  moment,  we  can  generalise  it  directly; 
therefore  in  the  rest  of  our  comments  we  treat  the  general  moment  rma  where  Besieris 
and  Tappert  treated  only  ru.  Their  equation  3.2  was  expressed  in  a  Fourier- 
transformed  domain,  but  can  be  expressed  in  our  notation  as 

«.rw.(x)  =  -.Lorw.(,)  -  /  (34) 

—  00 

This  equation  is  equivalent  to  (33)  to  order  L*M?,  and  it  should  be  noted  that  both  (33) 
and  (34)  are  invalid  if  LjM?  is  not  small.  Unlike  (33),  (34)  implies  a  “memory”  effect  in 
which  the  gradient  of  the  moment  depends  explicitly  on  the  moment  at  all  previous  «*s. 
The  Markov  approximation  leads  to  (7),  which  eliminates  the  memory  effect  and 
requires  only  a  correlation  function  of  the  medium  along  a  specified  (shifted)  direction. 
Besieris  and  Tappert  pointed  out  that  a  weaker  approximation,  called  the  “long-time 
Markov”  approximation  leads  to  a  local  (non-memory)  equation  (their  equation  3.3), 
that  in  our  notation  is  expressed  as 

dtrmm  -  -.t0r..(r)  -1/  < M[z) M(z') > |r„ (*) ,  (35) 


where  the  L0  operator  acts  only  on  Af(*f),  not  on  rM(r),  in  the  last  term.  We  have 
derived  (35)  by  use  of  the  Wigner-function  notation  of  Besieris  and  Tappert.  We  are 
only  considering  situations  in  which  the  parabolic  wave  equation  is  valid.  It  has  been 
shown  that  in  that  case  the  long-time  Markov  approximation  is  valid, M  and  therefore 
(35)  is  as  valid  as  (33). 

Because  L9  is  an  operator,  the  integrals  in  (33-35)  involve  the  medium  correlation 
function  in  all  directions,  or,  in  the  Fourier-transform  domain,  require  a  scattering  ker¬ 
nel  as  a  function  of  scattering  angle.  The  Markov  approximation  to  (33)  consists  of 
simplifying  the  deterministic  propagation  operator  e  for  *  -*'on  the  order  of 

Lr.  Instead  of  correlating  M[z)  with  all  possible  transverse  positions  of  Af(x*),  the 
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Markov  approximation  corresponds  to  choosing  only  one  transverse  position  for  Mi**). 
(See  Figure  1,  where  point  A  represents  an  arbitrary  transverse  position.)  If  the  wave 
represented  by  raa(s)  were  the  unperturbed  solution,  then  deterministic  propagation 
would  move  the  phase  in  the  direction  of  the  unperturbed  ray.  If  the  wave  energy  is 
travelling  close  to  the  unperturbed  ray  this  operator  retains  its  behavior  to  first 
approximation.  As  a  result,  deterministic  propagation  approximates  a  shift  along  the 
unperturbed  ray  to  point  B,  i.e.,  f*(**) «  1^,  ( t *),  where  the  ray  is  forced  to  go  through 
f*(r).  Hence  can  be  approximated  by  Mt^(ir).  This  is  the 

appropriate  definition  of  the  Markov  approximation  (rather  than  assuming  the  medium 
is  delta-correlated  along  the  t  axis)  and  it  immediately  yields  (7)  from  (33).  In  practice, 
instead  of  using  the  actual  unperturbed  ray,  the  tangent  to  the  ray  at  t  is  often  used. 

If  the  delta-correlated  assumption  were  made,  it  would  correspond  to  evaluating 
at  P0^  C,  which  is  strictly  valid  only  if  there  is  a  single  unperturbed  ray  trav¬ 
elling  along  the  i-axis.  If  the  medium  fluctuations  are  isotropic,  the  correlation  of  any 
point  at  *'  with  the  point  P  at  *  will  give  the  same  result  because  of  the  parabolic 
approximation,  and  hence  the  delta-correlated  assumption  is  as  good  as  any  other 
choice.  However,  for  an  anisotropic  medium  it  is  important  that  point  B  (and  hence 
(7))  be  used,  even  when  the  Markov  approximation  is  invoked.  Note  that  (7)  can  be 
used  in  the  presence  of  a  deterministic  background  refractive  index. 

The  difference  between  (33)  and  (7)  can  be  caused  by  directions  different  from  the 
unperturbed  ray  becoming  important.  It  is  in  this  sense  that  (33)-(35),  which  never 
refer  to  unperturbed  rays,  are  more  general  than  (7),  which  does.  A  transverse 
wavenumber  kT,  coming,  for  example,  from  M,  causes  the  angle  to  change  by 
Si  *  kT  /  k.  A  transverse  error  in  position  of  about  kTLP  /  k  is  made  by  assuming  the 
direction  of  the  unperturbed  ray.  Thus,  in  order  for  the  Markov  approximation  to  be 
valid,  it  is  required  that  kTLP  /  k  CLr,  where  LT  is  the  transverse  scale  of  concern. 
Since  LT  «  l/kT,  the  Markov  approximation  fails  at  sufficiently  small  k  «  LpfL\.  The 
parameter  a  »  kL$/Lp  introduced  by  Beran  and  McCoy^  and  discussed  further  in 
Flatte'  reflects  these  considerations.  For  small  a,  one  can  use  (33)  or  its  equivalent. 
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V.  Path-Integral  Derivation 

We  recapitulate  the  derivation  of  the  path-integral  expression  (7)  from  (10).  Using 
the  assumed  Gaussian  behavior  of  the  fluctuations,  we  obtain  from  (10) 

r„.  -  /  n  0 «?(*)•'  (36) 

where  50  is  the  part  of  5  in  (0)  that  does  not  involve  ft,  and 

V  -  -  *  **'  <  A f(*)A/(0  >  (37) 

m 

The  expression  (36)  is  an  exact  representation  of  the  moment  of  the  solution  of  the  par¬ 
abolic  equation  with  Gaussian  fluctuations.  It  is  not  used  in  practice  as  it  stands 
because  V  depends  on  the  paths  at  two  values  of  *,  namely  t  and 

The  Markov  approximation  for  the  path  integral  comes  from  assuming  that  the 
paths  do  not  stray  far  in  transverse  space  over  a  distance  lr\  they  all  move  approxi¬ 
mately  parallel  to  the  unperturbed  ray.  Thus,  in  the  Markov  approximation 

V  »  -jf  it  dz'  <  M{z)Mtiifl{z>)  >  (38) 

which  only  requires  knowledge  of  the  path  at  z.  The  final  result  (11)  follows  directly. 
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VL  Coordinate  Systems 

Moment  equations  can  be  formulated  in  a  variety  of  coordinate  systems,  while 
path  integrals  require  a  rectangular  coordinate  system.  There  has  been  a  fair  amount 
of  effort  expended  on  using  polar  coordinate  systems,  especially  for  point  source  prob* 
lems. 

The  same  results  (for  point  sources  among  others)  can  be  obtained  in  either  polar 
or  rectangular  coordinates.  Thus,  the  results  of  Shishovl14!  on  the  intensity  correlation, 
derived  in  spherical  polar  coordinates,  can  be  seen  to  be  identical  (after  an  appropriate 
transformation)  to  the  results  of  Codona  et  a/J16l  derived  in  rectangular  coordinates.  It 
was  necessary  for  Shishov  to  make  small  angle  approximations  in  addition  to  the  para* 
bolic  approximation  of  dropping  the  second  derivative  in  the  propagation  direction, 
whereas  Codona  et  al.  only  require  the  single  parabolic  approximation. 
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VJL  Summary 

We  have  derived  moment  equations  in  coordinate  reprenentation  under  the  Mar¬ 
kov  approximation  that  apply  in  anisotropic,  inhomogeneous  media  with  deterministic 
background.  The  derivation  shows  the  relationship  between  these  moment  equations 
and  modified  equations  that  are  valid  under  approximations  weaker  than  Markov;  the 
second-moment  equation  of  Besieris  and  Tappert  is  a  special  case  of  these  modified 
equations. 

In  a  hierarchy  of  approximations  we  begin  with  the  parabolic  wave  equation  itself. 
A  path  integral  with  non-local  exponent  can  be  written  as  an  exact  solution,  although  it 
is  not  yet  useful  in  practice.  The  next  level  is  the  approximation  that  the  interaction 
strength  over  a  correlation  length  is  small— this  “first-order  perturbation  theory”  leads 
to  the  modified  moment  equations,  and  in  homogeneous,  isotropic  media,  to  the  stan¬ 
dard  moment  equations  and  path-integral  expressions.  In  anisotropic,  inhomogeneous 
media,  however,  a  further  approximation  is  necessary  to  obtain  the  moment  equations 
and  path  integral  expressions.  This  further  approximation  is  that  the  significant  flow  of 
wave  energy,  or  the  important  paths,  are  parallel  to  the  unperturbed  ray;  we  call  this 
the  Markov  approximation  because  its  violation  implies  the  appearance  of  correlations 
between  successive  scatterings.  We  have  shown  that  the  moment  equations  and  the 
path-integral  expressions  for  the  moments  are  mathematically  equivalent  under  the 
Markov  approximation.  Thus  the  two  formalisms  have  exactly  the  same  physical  con¬ 
tent.  In  an  anisotropic  medium,  the  moment  equation  involves  a  shift  operation  to  cal¬ 
culate  the  medium  correlation  function  along  the  unperturbed  ray;  this  form  of  the 
moment  equation  has  not  been  given  before. 

We  have  also  pointed  out  that  all  appropriate  formulae  can  be  derived  in  a  rec¬ 
tangular  coordinate  system  (even  for  point  sources). 
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Appendfac 

We  must  show  that  the  scaling  pJj  |  ~  (Sx)%  bolds  for  integrals  of  the  form 

(a.i) 

If  F  b  expandable  in  a  power  series  (even  if  the  radius  of  convergence  is  zero)  this  result 
follows  immediately.  One  expands  F  and  integrates  term  by  terra,  obtaining  a  power 
series  in  ( Sx)\  By  standard  methods  in  the  theory  of  asymptotic  expansions,  only  the 
low  order  terms  need  to  be  retained  as  8z— *0. 

For  singular  functions,  a  demonstration  is  not  as  simple.  One  may  worry  about 
cancellations  between  terms  in  the  exponent,  since  the  signs  might  differ. 

We  will  content  ourselves  with  a  demonstration  in  the  case  likely  to  arise  in  prac¬ 
tice.  It  is  common  to  model  a  random  medium  as  having  a  power  law  structure  func¬ 
tion.  Thus  as  two  x’s  become  equal,  a  singularity  |z,  —  x,\r  with  p  >0  might  occur  in 
the  integrand.  In  order  to  have  possible  cancellations  in  the  exponent,  we  assume  that 
ki  **kj  =>k,  and  the  exponential  factor  is  exp (i*  [6z?  -  6xf)/26t ).  We  assume,  for  simpli¬ 
city,  that  *i  and  are  one-dimensional;  higher  dimensional  singularities  are  effectively 
weaker. 

Define  v  -  (Sx{  +  6x^/2,  p  —  -Say,  <*  —  ay  -sf.  The  singularity  from  the  pre¬ 
vious  step,  Xj  -  Xj  -  6xj  is  |p  -  op.  The  integral  to  be  evaluated  is 

/  in  iv  t<k,u’lh  |p  -  op  (A.2) 

We  would  like  to  ignore  the  p  dependence  in  f.  However,  spurious  large-p  contri¬ 
butions  would  arise,  even  though  we  are  only  interested  in  contributions  from  p  close  to 
a.  To  drop  the  p  dependence  of  f  and  also  to  simplify  the  analysis,  we  introduce  a  con¬ 
vergence  factor  exp(-*(p8  +  i/^/Sa1  “*)•  As  long  as  p,y  ~  6x*,  this  factor  does  not  change 
the  integral  as  Sa-»0  (we  are  assuming  c  >0).  Conversely,  if  the  integral  in  the  limit 
St— >0  does  not  depend  on  a  and  c,  then  p  and  v  are  of  order  6x*. 

The  integral  is  then 

I  m  J  ip  it/  e*'"'/1’  exp[  -a(p8  +*'*)/$*'  “*]  |p  -op  (A.3) 


The  p  integral  can  be  done: 


'  V»  '  ■■*.  ■"-*»  '*C®  *  *.  *.  ._  »  ^  •  «.  • 


I  m  c,  /  i¥  #••’/*•  «xp[  -  a(«r*  +  I^J/ll1  “*] 

,1  ,.l£±5iilU2£j 


(A.4) 


where  tfiii  confluent  hypergeometric  function  and  Cx  is  a  constant  independent  of 
fiyt,  and  a  (as  are  Ct  and  Ou  below).  The  hypergeometric  function  has  a  part  that 
behaves  as  the  exponential  of  its  argument  for  large  (positive)  values  of  its  argument,  a 
part  that  falls  as  a  power  (since  (?  + 1)/2  is  positive)  and  a  part  at  small  values  of  the 
argument.  These  last  two  parts  can  be  combined  into  a  bounded  part.  We  show  that 
the  exponential  part  gives  the  leading  behavior  and  the  bounded  part  is  a  higher  power 


of  St. 


The  contribution  /,  from  the  exponential  asymptotic  part  of  M  is 


It  -  C9  J  iv  eaar/fc  exp[  -  a(o*  +  rf/Sx'  -«] 


f(o, v)  a  +r)(1  “*!/*  exp  |- 


26  x*  >  or  -  it  t/V 
4a**1-* 


The  exponential  from  M  cancels  much  of  the  first  two  exponentials: 


l|  “  Ctjivtxp 

,  i 


!(a,*-)  +r)/»  Stf1  +*H1  -«Va  (A.6) 


which  can  be  done  explicitly.  Only  the  first  term  in  the  exponential  survives  as  St-*  0. 
The  result  is  independent  of  a  and  c,  and  is 

/,  -  St  fl«,0)M'  (A.7) 

exactly  as  would  be  obtained  from  the  Taylor  series  expansion  for  I. 

We  now  turn  to  the  contribution  1%  from  the  bounded  part  of  hi.  We  show  /s  has  a 
higher  power  of  St  than  /,.  We  can  set  <  to  any  positive  value.  At  large  <  we  depend  on 
the  fact  that  «<*"'’/*  averages  to  xero  for  v  ~  St1  +t  for  any  positive  6,  but  it  would  be 
necessary  to  examine  the  detailed  behavior  of  M  to  use  this  fact.  On  the  other  hand, 
for  small  enough  c,  it  suffices  to  bound  the  integral  by  the  integral  of  the  absolute  value 
of  the  integrand.  The  convergence  factor  provides  a  cutoff  at  v  ~  Sx1  ~‘/a.  Thus 
/  ivtxp\  -av*/!*1  “*]  f(o,i>)  gives  a  contribution  scaling  like  St1  -</2.  Thus  /8  is 
bounded  by  an  expression  which  scales  as 

...  .  _  ,  .  ..  ...  ... _  _ _ _ _ _ _ .VU-S5. 


/,  ~  -0/a  +(» +»)(*  -‘V*  .  ft( i  ♦>/*)(»  -«)  (A.8) 

As  long  ss  we  have  chosen  *  small  enough,  (he  exponent  of  ft  is  larger  than  1,  and  /* 
can  be  neglected  relative  to  /t.  Thus  we  have  established  the  necessary  scaling  of  p  and 
v  even  in  the  singular  case. 
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Figure  Captions 


Figure  la.  Moment-equation  expression  of  the  Markov  approximation.  The  corre¬ 
lation  should  be  taken  between  a  point  at  t  (point  P)  and  an  arbitrary 
point  at  s'  (point  A).  Instead  it  is  taken  with  the  point  B,  obtained  by 
extrapolating  along  the  unperturbed  ray  from  P.  The  assumption  of 
delta-correlated  medium  fluctuations  leads  to  the  incorrect  formulation 
of  correlations  between  points  P  and  C.  The  dashed  lines  indicate  the 
idea  of  a  scattering  as  a  function  of  angle  from  point  P. 


Figure  lb.  Path-integral  expression  of  the  Markov  approximation.  The  general 
path  at  *'  (point  A)  is  approximated  by  the  path  at  z  extrapolated 
along  the  unperturbed  ray  (point  B). 


Figure  2a.  Typical  z  values  of  the  interactions  from  a  Taylor  series  term  in  (9)  are 
indicated  by  x’s.  Dashed  lines  show  which  interactions  are  correlated. 
It  is  assumed  that  LfM?  Cl. 


Figure  2b.  A  portion  of  a  contribution  to  (9)  which  is  improperly  ordered  in  “first 
order  perturbation  theory.*’  Such  contributions  are  small  if  Cl. 
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A  relativistic  derivation  of  the  Schr&dinger  equation  is  given  from  principles 
known  to  physicists  in  1926.  Schrddinger's  preference  for  Hamilton's  optical* 
mechanical  analogy  over  the  relativistic  route  is  discussed  The  derivation  is 
given  of  a  classical  analog  to  the  Schr&dinger  equation  called  the  parabolic  wave 
equation,  which  describes  waves  propagating  in  a  narrow  angular  cone;  the  fact 
that  this  classical  version  was  not  discovered  until  about  1950  is  discussed 


Introduction 

Two  great  triumphs  of  twentieth  century  physics,  relativity  and  nonrela- 
tivistic  quantum  mechanics,  have  been  part  of  the  undergraduate  physics  00171* 
culum  for  only  a  few  decades.  Both  these  subjects  have  their  paradoxical 
aspects  with  which  students  must  grapple. 

It  is  generally  agreed  that  quantum  mechanics  is  more  difficult  to  accept 
than  relativity,  because  of  its  apparent  violations  of  cherished  ideas  such  as 
causality  or  locality.  However,  one  aspect  of  quantum  mechanics  that  should 
not  remain  mysterious  is  the  derivation  of  Schr5dinger*s  wave  equation.  On  a 
suggestion  of  Debye,  Schrodinger  set  himself  to  writing  a  wave  equation  for  an 
electron  around  a  proton.1  Once  given  the  idea  that  the  electron  might  be 
represented  by  a  wavefunction.  this  is  a  problem  in  classical  physics;  yet 
Schrddinger  came  out  with  a  "wave  equation”  that  no  one  had  seen  before. 

The  purposes  of  this  article  may  be  stated  as  follows:  1)  to  show  a  simple 
derivation  of  the  Schrodinger  equation  starting  from  a  classical  wave  equation 
and  some  physical  assumptions  that  would  have  been  plausible  to  physicists  in 
1926;  2)  to  discuss  why  Schrodinger  preferred  to  use  Hamilton's  optical* 
mechanical  analogy  rather  than  follow  the  relativistic  route  to  his  equation;  and 
3)  to  discuss  how  physicists  have  used  an  analog  to  the  Schrddinger  equation 
called  the  parabolic  wave  equation  in  completely  classical  contexts,  and  why 
nineteenth-century  physicists,  who  certainly  had  the  mathematical  tools,  did 
not  write  down  a  Schrddinger  equation  in  solving  some  classical  wave- 
propagation  problems.  Finally,  the  relationships  between  the  full  wave  equation, 
the  parabolic  wave  equation,  Huygens'  construction,  the  Schrddinger  equation, 
and  Feynman’s  path  integral  are  briefly  discussed. 


1.  The  Relativistic  Equation 

Electromagnetic  waves  obey  the  classical  wave  equation 

+r  -  £-«„/•  =  o  { t ) 

c 

where  F  may  be  a  component  of  the  electric  or  magnetic  field.  If  we  desire  the 
equation  that  can  represent  matter  as  a  wave,  following  de  Broglie,  we  must  deal 
with  the  problem  of  rest  mass.  Schrodinger  himself  solved  this  problem  in  1925, 
but  did  not  publish  until  1926  for  reasons  we  will  discuss  later.1  By  that  time 
Klein8  and  Gordon9  had  derived  the  same  equation  by  generalizing  Schr5dinger*s 
nonrelativistic  equation. 

The  matter-wave  equation  should  he  able  to  describe  something  like  a  parti¬ 
cle  at  rest  with  finite  energy  me2.  The  connection  between  energy  and  fre¬ 
quency  b>o  =  E/It  was  not  only  known  but  used  extensively  by  1926,  so  it  is  plau¬ 
sible  to  want  a  wave  function 

«  =  exp  -i  ~~t  (2) 

(with  no  space  dependence)  to  satisfy  the  fundamental  equation.  This  wave 
function  was  explicitly  suggested  by  de  Broglie  in  1925,  but  be  failed  in  his 
attempts  to  find  the  appropriate  equation.4  It  is  plausible  to  simply  add  a  covari¬ 
ant  term  to  (1)  in  order  to  have  (2)  satisfy  the  equation.  One  is  then  led  to  try 

V*u  -  jftu  u  -  =  0  (3) 

as  the  fundamental,  relativistic,  equation  for  free,  massive  particles.  This  is  the 
Klein-Gordon  equation.  The  wavefunction  u  has  something  to  do  with  a  free  par¬ 
ticle,  but  we  will  avoid  as  much  as  possible  discussing  interpretations  of  the 
wavefunction.  Equation  (3)  is  often  justified  as  the  wave  analog  of  the  relativistic 
energy-momentum  relation,  E 2  =  p8c8  +  m8c\  transformed  using  de  Broglie’s 
relations,  but  the  above  derivation  is  preferable  pedagogically.  The  addition  of 
the  new  m8  term  is  crucial  to  all  that  follows,  since  without  it  all  waves  would 
move  at  speed  c,  and  there  would  then  be  no  possibility  of  a  nonrelativistic 
approximation. 


0.  Derivation  of  Um  ftorPulicte  Schrftdtng «r  Equation 

We  desire  solutions  to  (3)  to  represent  particles  moving  through  space.  We 
have  one  solution  (Eq.  (2)),  but  it  is  not  too  interesting  since  it  describes  a  free 
particle  at  rest  whose  wavefunction  fills  all  of  space  and  oscillates  very  fast.  Let 
us  search  for  more  interesting  solutions  that  vary  in  both  space  and  time.  We 
try 


u  =  Hz,  t )  exp  I  - 


Putting  (4)  into  (3)  results  in 


+*-${*  -  i-r ♦  ^ = °  •  <*> 

This  is  a  general  equation  for  if.  Since  it  is  an  equation  that  is  linear  in  f,  it 
makes  sense  to  consider  each  Fourier  component  of  if  separately: 

i  ~  exp  |i(fc*  -  ut)  J  (6) 

where,  of  course,  k  =  2n/A  and  A  is  the  wavelength.  If  we  restrict  ourselves  to 
low-frequency  solutions  in  which 

»«  (7) 

then  the  d«  term  in  (5)  is  negligible,  and  (5)  becomes 

(8) 

which  is  the  free-particle  Schr&dinger  equation.  Note  that  the  first  derivative  in 
time  and  its  imaginary  coefficient  come  naturally  from  the  second  derivative 
because  the  time  dependence  of  if  is  a  difference  frequency  from  the  fundamen¬ 


tal  CJg. 


The  low-frequency  requirement  (7)  can  be  stated  in  two  ways.  The  first  is 


obvious: 


htt  «  me*  ; 


that  it,  interpreting  the  frequency  of  ^  as  an  energy,  that  energy  must  be  much 
less  that  the  rest-mass  energy.  A  second  interpretation  is  afforded  if  we  imagine 
making  a  wave  packet  of  waves  like  (6).  Equation  (8)  implies 


(10) 

(11) 

(12) 


In  other  words,  the  Schrodinger  equation  (8)  is  a  nonrelativistic  approximation 
to  the  Klein-Gordon  equation  (3).  Finally,  (9)  and  (10)  can  be  combined  to  give 
another  version  of  the  validity  requirement: 


A  » 


A 

me 


(13) 


In  other  words,  the  Schrddinger  equation  is  a  large-wavelength  approximation  to 
the  full  wave  equation;  that  is,  the  wavelength  must  be  large  compared  to  the 
Compton  wavelength  of  the  electron. 


m.  Addition  of  an  External  Potential 

The  addition  of  electromagnetic  potentials  to  (3)  should  follow  the  rules  of 
relativistic  invariance.  It  was  known  in  1928  that  f  and  £/c.  the  scalar  and  vec¬ 
tor  potentials,  are  components  of  a  four-vector.  The  minimal  coupling  implies 
that  (3)  becomes 

(v  -  Sr*'  -  ?-<»• -  r*)*  •  ir  M  = 0  <14> 

Consider  the  case  in  which  4  =  0.  In  that  case  the  equation  (5)  for  becomes 

**  -  ♦  »«  -  f-i-a.  -  =  o  (is) 
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The  size  of  atoms  (known  in  1928)  implies  that  the  relevant  Coulomb  potentials 
felt  by  an  electron  are  small  compared  with  me*;  that  is 

tip  «  me9  .  (16) 

Therefore  the  terms  in  (14)  satisfy  a  hierarchy  such  that  the  last  three  terms 
are  negligible.  Again  this  corresponds  to  i>  having  low  frequency;  that  is  (9)  and 
(12)  are  satisfied.  The  resulting  equation  is 

ihdtf  =  -  +  e  p  .  (17) 

that  is,  the  full  nonrelativistic  Schrodinger  equation  with  external  potential. 

IV.  The  Fine  Structure  of  Hydrogen 

Schrodinger  wrote  down  equation  (14)  in  1925,  but  did  not  mention  it  in  a 
publication  until  late  1926.  He  gave  as  his  reason  that  an  exact  solution  for  the 
energy  levels  of  hydrogen  disagrees  with  experiment  in  the  fine  structure.1 
Since  fine-structure  energy-level  differences  are  quite  small,  this  amounts  to 
saying  that  the  neglected  three  terms  in  (15)  are  incorrect.  The  actual  numbers 
are  elegantly  worked  out  in  ScbifTs  text,0  where  it  is  shown  that  the  hydrogen 
fine  structure  from  the  last  three  terms  on  the  left  of  (15)  is  about  twice  as 
large  as  experiment,  and  has  a  slightly  different  dependence  on  the  orbital 
quantum  number. 

Schr&dinger  himself  pointed  out  that  the  discrepancy  probably  bad  some¬ 
thing  to  do  with  Gouda mit  and  Uhlenbeck’s  hypothesis  of  electron  spin.0  This  is 
probably  the  main  reason  that  Schrodinger  did  not  derive  his  nonrelativistic 
time-dependent  equation  along  the  lines  of  our  discussion.  And  of  course  it  is 
one  of  the  reasons  why  Dirac's  contribution,  which  seemed  to  find  the  existence 
of  electron  spin  in  a  "natural”  way,  is  so  admired. 

However,  as  will  be  discussed  later,  the  derivation  of  the  Schrodinger  equa¬ 
tion  from  a  classical  wave  equation  is  much  to  be  desired  pedagogically.  Given 
the  many  paradoxes  of  quantum-mechanical  matter  waves  in  1926,  including  the 
existence  of  spin,  the  interpretation  of  the  wavefuncUon,  and  the  inability  to 
create  a  dispersion-free  wave-packet,  the  hydrogen-flne-structure  discrepancy 
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was  hardly  a  reason  to  abandon  the  whole  approach. 
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V.  Antimatter  and  Electron  Spin 

Dirac’s  prediction  of  antimatter  from  his  relativistic  equation  for  the  elec* 
tron  with  spin  is  rightly  considered  one  of  the  triumphs  of  modern  theoretical 
physics.  It  is  possible  to  make  the  same  point  from  the  Klein-Gordon  equation, 
by  simply  pointing  out  that  the  nature  of  equation  (3)-second  order  in  time, 
with  m*  appearing  rather  than  m ,  allows 

G  =  exp[  +  i^-/]  (10) 

to  be  a  solution  as  well  as  u  from  (2).  The  same  derivation  of  the  nonrelativistic 
equation  goes  through,  and  as  long  as  e<p  «  me2,  the  two  states  will  not  mix.  As 
the  energy  etp  gets  larger,  the  coupling  between  the  two  states  becomes  more 
important,  leading  in  lowest  order  to  vacuum  polarization,  Zitterbewegung,  etc., 
and  eventually  to  antimatter  production  and  the  necessity  for  second  quantiza¬ 
tion  and  quantum  electrodynamics. 

It  was  the  unwelcome  existence  of  negative-energy  solutions  to  the  Klein- 
Gordon  equation  that  led  Dirac  to  search  for  an  equation  in  first  derivatives.  His 
search  led  eventually  to  four  component  spinors,  and  thus  he  was  forced  back  to 
the  situation  that  is  evident  in  the  Klein-Gordon  equation  at  the  start— the  simul¬ 
taneous  existence  of  positive  and  negative-energy  solutions. 

The  discrepancy  in  the  hydrogen  fine  structure  cannot  be  dismissed  as  an 
effect  of  relativity.  The  fact  is  that  the  spin  of  the  electron  is  1/  2,  not  0.  Some 
people  feel  that  starting  with  the  Klein-Gordon  equation  to  describe  the  electron 
is  therefore  too  misleading.  It  should  be  pointed  out  that  learning  the  Klein- 
Gordon  route  is  quite  valuable,  since,  as  Case7  has  shown,  the  Foldy-Woutbuysen 
transformation  makes  the  resulting  forms  for  operators  and  the  nonrelativistic 
limits  quite  similar  for  integer  and  half-integer  spins. 
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VI.  Relation  to  Hamilton's  Optical-Mechanical  Analogy 

It  is  of  the  greatest  importance  to  realize  that  (17)  is  a  low-frequency 
approximation  to  the  relativistic  equation.  It  was  Schrddinger's  great  desire  to 
make  the  connection  between  quantum  mechanics  and  classical  mechanics  by 
the  Hamiltonian  method  that  Sommerfeld  and  Runge  used  to  connect  wave 
optics  to  ray  optics.8  That  connection  requires  taking  the  high-frequency  limit  of 
the  wave  theory.  Schrodinger  was  fully  aware  that  the  frequency  was  associated 
with  the  energy,  and  that  the  logical  total  energy  would  involve  a  much  larger 
value,  me2.  It  was  no  problem  for  him  to  realize  that  the  requirement  for  the 
classical  equations  of  motion  to  be  valid  would  be  a  requirement  that  the 
wavelength  is  small  compared  with  some  characteristic  length  L  for  the  poten¬ 
tial  to  change: 

X  «  L  .  (19) 

which  combines  with  (10)  to  give 

£-1?  «  »  .  (20) 


Thus  Schrodinger  took  the  point  of  view  that  classical  mechanics  should  be 
included  in  his  equation,  and  took  the  high-frequency,  small-wavelength  limit  to 
satisfy  that  requirement,  but  he  ignored  the  fact  that  the  nonrelativistic 
requirement  at  the  same  time  restricted  the  frequency  to  be  small  compared 
with  a  fixed  quantity  (me2/ ft).  He  did  realize  implicitly  that  a  restricted  range 
of  frequency  does  exist  that  satisfies  both  the  nonrelativistic  and  the  classical 
mechanics  requirements: 


ft 

2mLe 


«  u  « 


me2 

ft 


(21) 


Note  that  as  ft  -» 0  the  range  of  validity  of  nonrelativistic  classical  mechanics 
gets  larger  at  both  ends.  That  is  why  the  ft  -»  0  limit  is  often  used  as  the  classical 
limit. 

It  is  also  useful  to  put  the  requirement  (21)  in  terms  of  a  spatial  wavelength 
requirement  (X  =  2n/k); 


(22) 


The  left-hand  inequality  expresses  the  nonrelativistic  requirement  and  the  vali¬ 
dity  of  the  SchrSdinger  wave  equation,  while  the  right-hand  inequality  expresses 
the  classical-mechanics  requirement. 

Schrodinger’s  novel  use  of  Hamilton’s  '-'jthods,  wherein  he  used 
de  Broglie’s  relations  to  define  u  and  k  in  terms  of  E  and  p  and  then  created  a 
wave  equation  by  identifying  u  and  k  with  operators  on  a  wavefunction,  allowed 
him  to  bypass  the  requirements  of  relativity  and  allowed  him  to  ignore  the  left- 
hand  inequality  in  (22).  By  this  method  of  starting  with  classical  nonrelativistic 
equations  he  avoided  both  the  pitfalls  and  opportunities  associated  with 
antimatter  and  electron  spin. 


VII.  The  Classical  Schrodinger  Equation:  The  Parabolic  Wave  Equation 

The  Schrodinger  equation  comes  out  of  the  Klein-Gordon  equation  solely 
because  of  the  m2  term.  If  the  m2  term  were  not  there,  all  waves  would  move  at 
speed  c ,  and  there  would  be  no  hope  for  solutions  with  small  velocity.  Therefore 
the  classical  wave  equation  (1)  can  never  yield  the  Schrodinger  equation  with  a 
first  derivative  with  respect  to  time. 

However,  an  analog  of  the  Schrodinger  equation  in  which  the  first  derivative 
is  with  respect  to  a  spatial  coordinate  can  be  obtained  in  the  following  way.  Sup¬ 
pose  the  speed  of  the  waves,  c ,  is  a  function  of  position,  and  we  look  for  time- 
harmonic  solutions 

F(*,t)  =  u(*)e^u  (23) 

Then  our  equation  becomes  the  Helmholtz  equation 

V*u  +  k*u  =0  (24) 

where  k  =  u/  c  is  the  wavenumber,  which  is  a  function  of  position  because  of  c . 

The  general  solution  to  this  equation  involves  waves  in  all  directions.  How¬ 


(25) 


ever,  let  us  single  out  the  z  direction  as  of  special  interest,  and  write 

[V?  +  +  Jfcz]ii  =  0 
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where  Vf  =  dm  *  Ow.  Equation  (25)  hat  tome  of  the  characteriatica  of  the 
Klein-Cordon  equation  (3).  Let  ua  try  a  aotution  of  a  form  which  reminda  us  of 
waves  travelling  in  the  positive  s  direction: 

«  =  (28) 

The  resulting  equation  for  f(z,  y,  a)  is 

Vfy  +  [8*.  +  2 i*08.  +  =  0  (27) 

which  is  remarkably  analogous  to  (5).  We  see  that  the  term  is  playing  the 
role  of  the  m*  term  in  the  Klein-Gordon  equation.  If  kz  is  slowly  varying  in 
space,  then  we  can  pick  k0  such  that 

JJfcf  -kz\  «  kl  (28) 

In  other  words,  the  wave  speed  c  (g)  varies  by  a  small  fraction  of  itself,  and  does 
so  slowly  with  respect  to  g.  In  that  case,  i*  will  have  only  small-wavenumber 
components 

1 Ks)  ~  «<X  *  :  |q|  «  *0  (29) 

and  the  0„  term  in  (27)  is  negligible.  The  resulting  equation  is 

2tM,f  =  -VfV'  +  {let  -kz)f  ,  (30) 

which  is  exactly  analogous  to  the  SchrSdinger  equation  (17),  except  that  V?  is  a 
two-dimensional  Laplacian  rather  than  a  three-dimensional  one.  The  constants 
in  (30)  are  a  bit  different  from  (17).  The  analogy  more  direct  Let 

*(*)==  *o[l  +  M)  (31) 

where  /4(g)  is  the  variation  from  unity  of  the  index  of  refraction.  Now  write  (30) 
as 


and  write  (17)  as 


2  =  -Vft  +  2  kf/tf 


(32) 


-■  -v-.v-v  . 
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(33) 

The  analogy  is  complete,  and  we  see  that  the  appropriate  £q  is  the  inverse 
of  the  Compton  wavelength  of  the  electron,  and  the  variable  index  of  refraction 
is  analogous  to  the  potential  V  as  a  fraction  of  me*. 

There  is  an  analogy  to  antimatter  in  (32).  We  can  find  solution*  to  (25) 
which  are  travelling  in  the  negative  z  direction.  These  cure  analogous  to  .  ntipar- 
ticle  solutions.  As  long  as  p  is  slowly  varying  and  small,  the  waves  going  in  the 
positive  and  negative  directions  are  decoupled.  However,  if  fi  is  strong,  or  varies 
over  a  small  distance,  then  coupling  will  occur  in  the  form  of  backscattering. 
This  is  a  classical  realization  of  Feynman’s  picture  of  positrons  being  electrons 
moving  backwards  in  time.® 

Equation  (32)  is  called  the  parabolic  wave  equation  in  the  classical  context 
in  which  it  is  used.10  It  is  most  used  in  solving  problems  in  wave  propagation 
through  continuously  variable  media;  for  example  light  through  the  atmosphere, 
radio  waves  through  the  ionosphere  or  interplanetary  plasma,  or  sound  through 
the  ocean.11 

Let  us  discuss  the  validity  requirements  for  the  parabolic  wave  equation. 
The  main  ingredient  is  (28)  which  states  that  the  variations  in  the  wave  speed 
must  be  small  compared  with  unity,  in  analogy  with  the  nonrelativistic  require¬ 
ment  (16). 

In  the  quantum-mechanical  context  the  nonrelativistic  requirement  can  be 
expressed  in  a  variety  of  ways;  for  example,  that  the  wavelength  be  large  com¬ 
pared  to  some  quantity.  The  analog  of  this  requirement  in  the  classical  context 
is  quite  different,  because  the  analog  of  the  wavelength  of  the  quantum- 
mechanical  is  not  the  classical  wavelength  of  the  propagating  wave. 

In  the  classical  context  we  must  distinguish  between  the  longitudinal  and 
transverse  components  of  q,  the  wavevector  of  i>,  and  furthermore  we  have  the 

M 

wavelength  X,  defined  from  Aq. 

Let  9  =  (kj-.ki)  where  kf  has  components  in  the  x-y  plane,  and  kL  is  along 
the  z  direction.  The  first  term  of  (32)  must  be  larger  than  9„f.  This  results  in 
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the  requirement 


*r  «  k o  (34) 

which  means  that  the  total  wavevector  &  is  directed  at  a  small  angle  to  the  z 

axis.  This  latter  requirement  is  the  analog  of  the  nonrelativistic  requirement, 
v  «  c  in  the  Schr&dinger  equation  case. 

The  relation  of  the  various  wavenumbers  to  the  characteristics  of  the 
medium  index  of  refraction  is  quite  involved,  and  lies  beyond  the  scope  of  this 
paper,  because  it  brings  in  the  strength  of  the  fluctuations  (ji)  as  well  as  the 
scale  L. 

The  analogy  between  the  parabolic  wave  equation  and  the  Schrodinger 
equation  is  expressed  pictorially  in  Figures  1  and  2,  and  is  summarized  in 
Table  1. 

Vm.  History  of  the  Parabolic  Wave  Equation 

The  parabolic  wave  equation  cc'ild  be  called  the  classical  Schrodinger  equa¬ 
tion.  The  strongest  difference  is  just  one  of  variables,  since  the  parabolic  equa¬ 
tion  is  first  order  in  a  space  derivative  rather  than  the  time  derivative.  A  para¬ 
bolic  wave  equation  was  written  down  for  the  first  time  in  1948  by  Leontovicb  and 
Fock. 19  but  they  did  not  really  write  down  (32).  They  wrote  down  a  parabolic 
equation  for  radio  wave  propagation  over  the  surface  of  the  earth,  with  no  poten¬ 
tial  term.  Their  effect  was  controlled  by  diffraction  due  to  the  curved  boundary 
condition  over  the  surface  of  the  spherical  earth. 

Fock  was  one  of  the  pioneers  in  searching  for  wave  equations  for  quantum 
mechanics,  so  it  is  not  surprising  that  he  was  the  first  to  point  out  the  relation¬ 
ship  between  the  parabolic  wave  equation  and  the  Schrodinger  equation.19  He 
also  pointed  out  that  the  dependence  of  the  index  of  refraction  for  radio  waves 
on  height,  which  occurs  because  the  atmosphere  is  stratified,  is  an  analog  to  the 
quantum  mechanical  potential.  It  is  perhaps  surprising  that  be  did  not  make  the 
connection  until  1950;  he  also  did  not  discuss  the  requirements  for  validity  as  is 
done  in  this  paper. 
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Thanks  to  Fock,  the  parabolic  wave  equation  with  external  potential  at  in 
(32)  was  known  to  Soviet  worker*  in  the  1960*s.  It  was  common  knowledge  by 
the  time  of  Tatarskii's  classic  monograph  on  light  through  the  turbulent  atmo¬ 
sphere.14  (Tatarskii  gives  neither  references  nor  a  very  straightforward  deriva¬ 
tion.)  A  further  interesting  case  in  which  the  potential  term  is  replaced  by  a 
nonlinear  term  was  written  down  by  Kelley,18  who  was  dealing  with  self-trapped 
laser  beams. 

One  may  also  ask  why  the  parabolic  wave  equation  was  not  written  down  one 
hundred  years  earlier,  since  all  the  mathematical  tools  were  available  once 
Hamilton  had  used  complex  numbers  to  represent  oscillatory  phenomena.  Why 
did  not  those  mathematical  physicists  interested  in  the  wave  theory  of  light 
write  down  a  parabolic  equation  for  the  propagation  of  light? 

The  parabolic  wave  equation  is  fundamentally  a  small-wavenumber  approxi¬ 
mation.  Most  physicists  were  working  with  the  large-wavenumber  limit  of  wave 
propagation,  namely,  geometrical  optics.  This  was  especially  true  of  Hamilton, 
whose  work  had  tremendous  influence  over  those  who  followed. 

Probably  more  importantly,  experiments  were  not  done  in  continuous 
media.  Typical  problems  involved  diffraction  around  obstacles,  and  for  these 
problems,  solutions  in  the  form  of  integrals  were  directly  written  down. 

Rayleigh  typifies  the  attitude  of  several  generations  of  physicists  by  refer¬ 
ring  to  the  problem  of  continuous  media  in  his  Theory  of  Sound;14 

The  variation  is  supposed  to  be  so  slow  that  no  sensible  reflection 
occurs,  and  this  is  not  inconsistent  with  decided  refraction  of  the  rays 
in  travelling  distances  which  include  a  very  great  number  of 
wavelengths...  The  further  development  of  this  part  of  the  subject 
would  lead  us  too  far  into  the  domain  of  geometrical  optics.  The  funda¬ 
mental  assumption  of  the  smallness  of  the  wavelength...,  having  a  far 
wider  application  to  the  phenomena  of  light  than  to  those  of  sound,  the 
task  of  developing  its  consequences  may  properly  be  left  to  the  cultiva¬ 
tors  of  the  sister  science. 
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in  one  paragraph  be  washes  his  hands  of  the  question  and  tells  his  optical 
colleagues  that  it  is  a  geometrical  optics  problem)  It  was  not  until  1948,  when 
Feynman17  wrote  down  his  path  integral,  that  the  relation  between  the  various 
integral  solutions  to  the  differential  equations  describing  wave  propagation  were 
made  clear. 

DC  Huygens-ftresnel-Feynman  Theory 

Huygens  attempted  to  describe  light  as  a  wave,  in  analogy  to  sound.  Ironi¬ 
cally,  he  used  his  ideas  about  secondary  wavelets  to  prove  that  light  travelled  in 
straight  lines;  that  is,  he  dealt  only  with  geometrical  optics.1*  In  fact,  Huygens* 
construction  is  mathematically  equivalent  to  the  small-wavelength  limit  of  the 
full  wave  equation,  namely  the  unfolding  of  a  contact  transformation  (the 
Hamilton-Jacobi  equation  for  the  phase).1* 

Over  one  hundred  years  passed  before  attempts  were  made  to  generalize 
Huygens  construction  to  include  diffraction.  Fresnel’s  first  attempt  to  do  so 
uncovered  disturbing  difficulties  such  as  the  backward  wave  and  the  require¬ 
ment  that  the  secondary  wavelets  had  to  be  emitted  one-quarter  cycle  out  of 
phase  with  the  incident  wave.80  Helmholtz  and  Kirchhoff  attempted  to  improve 
on  Fresnel  by  introducing  the  so-called  obliquity  factor  (which,  among  other 
things  eliminated  the  backward  wave),  but  neither  they  nor  anyone  else  have 
been  able  to  create  a  Huygens-like  construction  that  provides  a  solution  to  the 
full  wave  equation.21 

After  Fresnel,  over  one  hundred  more  years  passed  before  Huygens’  con¬ 
struction  was  generalized  in  a  way  that  exactly  solved  a  wave  equation. 
Feynman’s  path  integral  is  an  exact  solution  of  the  Schro  dinger  equation.82 
Therefore,  when  translated  into  classical  terms,  it  is  an  exact  solution  of  the 
parabolic  wave  equation.  No  obliquity  factor  appears,  and  the  quarter-cycle 
phase  shift  is  relegated  to  the  status  of  a  mathematical  normalization  factor. 
No  backward  wave  appears,  since  no  coupling  to  a  backward  wave  is  allowed  in 
the  parabolic  wave  equation,  as  discussed  above.  Thus  the  Feynman  path 
integral  can  be  regarded  as  the  simplest  Huygens-Fresnel  theory  repeatedly 
applied  over  infinitesimal  steps  in  range. 
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Feynman  «u  concerned  with  quantum  mechanic*,  not  classical  wave  propa¬ 
gation;  hence  he  worked  with  the  Schrbdlnger  equation.  The  analog  of  generalis¬ 
ing  to  the  full  classical  wave  equation  would  have  been  to  try  to  solve  all  the 
problems  associated  with  relativity,  antimatter,  and  particle  production.  Thus 
the  road  of  Helmholtz  and  Kirchhoff  was  not  particularly  tempting  in  the 
quantum-mechanical  context.  One  could  say  that  Fresnel's  difficulties  were 
solved  not  by  constructing  a  solution  to  the  full  wave  equation,  but  by  going  back 
to  Fresnel’s  simple  construction  and  realizing  that  it  was  the  solution  to  an 
approximate  wave  equation,  one  that  is  valid  for  waves  propagating  in  a  narrow 
angular  cone. 

Geometrical  optics  is  valid  in  the  limit  of  infinitesimal  wavelength.  The  par¬ 
abolic  wave  equation  and  the  Feynman  path  integral  are  valid  for  waves  pro¬ 
pagating  in  a  narrow  angular  cone. 

X  Conclusion 

In  1926,  the  Schr6dinger  equation  was  a  wave  equation  that  had  never  been 
seen  before.  Yet  it  does  appear  in  classical  contexts,  a  fact  only  realized  well 
after  Schrodinger's  work.  In  classical  form  it  is  called  the  parabolic  wave  equa¬ 
tion,  and  its  derivation  from  the  full  classical  wave  equation  involves  a  restric¬ 
tion  to  waves  travelling  within  a  narrow  cone  of  angles  in  a  particular  direction. 
An  analogous  derivation  of  SchrSdinger's  equation  begins  with  the  Klein-Gordon 
equation,  and  need  only  require  that  the  potential  energy  is  small  compared 
with  me8.  Undergraduates  might  well  be  introduced  to  Schrodinger's  equation 
by  this  route  rather  than  the  standard  ad  hoc  approaches. 

The  development  of  theories  of  wave  propagation  through  continuous 
media,  sometimes  with  a  random  component,  owes  much  to  the  use  of  the  para¬ 
bolic  wave  equation.  The  analogy  between  that  equation  and  the  Schrddinger 
equation  has  brought  the  full  power  of  much  that  was  learned  in  nonrelativistic 
quantum  mechanics  to  bear  on  classical  wave  propagation  problems.1*  Most  not¬ 
ably,  Feynman's  path-integral  technique  has  been  important  in  solving  problems 
in  wave  propagation  in  random  media  (WPRM).  Perhaps  in  the  future,  quantum- 
mechanical  problems  in  condensed  matter  may  be  helped  by  recent  progress  in 
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WPRM,  where  the  medium  to  considered  as  having  a  statistically  fluctuating  wave 
speed.  That  would  have  been  appreciated  by  Schrbdinger. 
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Figure  Captions 


1.  Typical  solutions  to  classical  wave  equations. 

(a)  The  full  wave  equation  allows  waves  to  travel  in  all  directions.  For  an 
incident  plane  wave,  such  solutions  will  be  generated  by  medium  variations 
with  scales  comparable  with  a  wavelength.  Waves  at  large  angles  to  the 
forward  or  backward  directions  can  be  thought  of  as  linear  combinations 
(couplings)  of  nearly  forward  and  backward  waves. 

(b)  The  parabolic  wave  equation  describes  only  waves  travelling  in  the 
nearly  forward  or  backward  directions,  and  allows  no  coupling  between 
them. 

2.  Typical  solutions  to  quantum-mechanical  wave  equations. 

(a)  The  Klein-Gordon  equation  allows  particle  and  antiparticle  waves  to 
travel  at  speeds  up  to  the  speed  of  light.  For  an  incident  plane  wave,  such 
solutions  will  be  generated  by  external  potentials  that  have  strengths  com¬ 
parable  with  me*.  Waves  at  high  speeds  can  be  thought  of  as  linear  combi¬ 
nations  of  particle  and  antiparticte  waves,  representing  the  possibility  of 
(real  or  virtual)  particle  production. 

(b)  The  Schr&dinger  equation  describes  only  slowly  moving  particles  or 
antiparticles,  and  allows  no  particle  production. 
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a)  FULL  WAVE  EQUATION 


b)  PARABOLIC  WAVE  EQUATION 
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I  provide  here  e  guide  to  a  few  recent  result#  obtained  by  the  group  at  the  Center 
for  Studies  of  Nonlinear  Dynamics  in  the  area  of  wave  propagation  in  random  media. 


I.  Average  Arrival  Time  of  Wave  Pubes  Through  Continuous  Random  Media 


A  random  medium  consisting  of  discrete  point  scatterers  in  a  homogenous  back¬ 
ground  will  delay  the  arrival  of  a  pulse,  because  all  scattered  paths  are  of  greater 
length  than  the  straight-line  path  from  transmitter  to  reciever.  A  continuous  random 
medium  is  different. 

Our  recent  analysis1  points  out  that  a  fluctuating  continuous  medium  can  cause  an 
average  advance  of  the  pulse  arrival  time.  All  previous  analyses  have  dealt  with  situa¬ 
tions  in  which  pulses  are  delayed  on  the  average.2,3  By  convention,  the  ensemble  aver¬ 
age  of  a  random  medium  is  taken  as  the  medium  reference  state,  and  the  small  fluctua¬ 
tions  about  this  reference  state  are  thus  by  definition  a  zero-mean  1  andom  process.  The 
arrival-time  advance  or  delay  is  relative  to  the  travel  time  through  the  reference  state. 
Thus,  for  example,  results  through  turbulent  air  or  plasma  are  relative  to  quiescent  air 
or  plasma,  not  vacuum. 

The  behavior  of  a  wave  propagating  through  a  random  medium  is  controlled  by 
relationships  between  the  wavenumber  (k)  of  the  propagating  wave,  the  range  (If),  and 
the  strength  and  size  of  the  medium  fluctuations.4  Unsaturated  behavior  corresponds  to 
one  stationary-phase  path  (ray),  and  occurs  if  the  medium  fluctuations  are  weak 
enough.  In  fully  saturated  behavior  the  original  ray  breaks  up  into  many  new  microrays 
which  are  statistically  independent  of  each  other.  Propagation  through  a  medium  of 
discrete  scatterers  falls  in  this  category.  Partially  saturated  behavior  occurs  in  a 
strongly  fluctuating  medium  with  a  power-law  spectrum,  which  has  enough  small-scale 
fluctuations  to  cause  the  breakup  into  many  microrays,  and  enough  large-scale  fluctua¬ 
tions  to  make  the  microray  bundle  behave  like  a  single  ray  in  its  wandering  from  the 
unperturbed  ray.  Experiments  in  waves  propagating  through  continuous  random  media 
typically  fall  into  this  category.  We  deal  only  with  the  important  case  in  which  the 
transverse  wandering  from  the  unperturbed  ray  is  small  compared  with  the  range  of 
propagation. 
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Briefly  our  results  are  as  follows:  if  the  travel  time  of  a  pulse  is  averaged  over  an 
ensemble  of  the  random  medium,  with  each  pulse  weighted  by  its  intensity,  then  the 
average  pulse  is  delayed,  regardless  of  the  type  of  propagation  behavior,  in  agreement 
with  previous  results.3,3  However,  if  the  average  travel  time  is  obtained  without  weight¬ 
ing  by  pulse  intensity,  then  a  pulse  advance  is  expected  for  both  unsaturated  and  par¬ 
tially  saturated  behavior,  while  a  pulse  delay  remains  for  the  fully  saturated  case.  The 
difference  between  intensity-weighted  and  unweighted  travel  time  probes  the  variance 
of  the  first  derivative  of  the  refractive  index,  smoothed  over  a  microray  bundle. 

The  effect  is  illustrated  by  a  simple  special  case.  Consider  a  point  source  and 
point  receiver  separated  by  range  R,  and  a  homogeneous  medium  in  the  absence  of 
fluctuations,  so  that  the  unperturbed  ray  from  source  to  receiver  b  a  straight  line.  The 
random  medium  is  concentrated  in  a  “phase  screen’*  at  a  distance  z  from  the  source. 
Thb  screen  has  the  effect  of  advancing  the  time  of  a  wavefront  by  a  random  amount 
t(x)  where  x  is  the  position  on  the  screen,  and  ((i)  b  a  stationary  Gaussian  random 
process  with  zero  mean.  (We  take  x  as  one-dimensional  for  simplicity.) 

Weak  fluctuations— In  the  geometrical  optics  limit  only  one  ray  exbts  from  source  to 
receiver.  The  travel  time  for  a  path  through  point  x  b: 

T{x )  «  +  O.Seol  Ax 3  -  t{x)  (1.1) 

where  A~l  =  z(R  —  z)/R.  By  Fermat’s  principle  the  ray  b  at  the  point  xt  such  that 
r(zr)  b  a  minimum.  For  the  case  of  weak  fluctuations  we  may  expand  t(x)  as 

t{x)  -  40  +  l'* +0.54"*3  (1.2) 

The  position  of  the  ray  follows  to  first  order  as 

xt  =  A~lc0t'  (1.3) 

The  travel  time  of  the  ray  b  then 

T(xr)  =  r0  +  O.5c0^-I<'3-/0-c0^-,<'3  (1.4) 

Thb  case  requires  that,  typically 

M-1**!  <  ig  (i.5) 

But  t  (and  hence  t0  and  4*)  are  (by  construction)  random  variables  with  zero  mean. 


Therefore  (he  t9  term  will  disappear  in  (he  average  (ratel  (ime  and  (he  only  effee(  of 
(he  fluctuations  will  come  from  (he  t*  (erma.  These  terms  arise  because  (he  ray  has 
moted  away  from  its  unperturbed  position.  The  first  V  term  is  positive,  corresponding 
to  a  pulse  delay,  and  represents  (he  effects  of  geometry;  the  perturbed  path  is  physically 
longer  than  (he  unperturbed  one.  The  second  t'  term  is  negative,  corresponding  to  a 
pulse  advance;  we  call  this  the  Fermat  term;  the  ray  sought  out  a  region  of  the  medium 
with  a  pulse  advance.  The  Fermat  term  is  twice  as  large  in  magnitude  as  the  geometry 
term.  The  average  travel  time  is 

<r>  -  r0  —  o.5c0a-1  <trt>  (1.6) 

so  that  the  pulse  on  the  average  arrives  early. 

There  is  a  subtlety  to  this  result.  In  the  we&k-fiuctuation  limit  the  intensity  is 
controlled  by  the  focussing  due  to  the  curvature  of  the  wavefront  as  it  exits  the  phase 
screen.  It  is  not  difficult  to  show  that  the  intensity  /  is,  to  first  order, 

/  =  l+X-lc0«»  (1.7) 

Consider  the  intensity-weighted  average  travel  time: 

<  JT(xf)  >  -  T0  -  0.5 c0A~l  <tn>  -  c9A~ 1  <l0t«>  (1.8) 

where  the  last  term  comes  from  the  correlation  between  the  intensity  and  the  travel 
time. 

For  any  random  function  l(x)  whose  Fourier  components  are  uncorrelated  (i.e.,  the 
correlation  function  is  translation-invariant): 

<t9t">  =  -<!**>  (1.9) 

Therefore 

<  IT(xr)  >  *  T0+0.5c0A-|<lrt>  (1.10) 

In  other  words,  the  intensity-weighted  average  travel-time  is  delayed  by  fluctua¬ 
tions  by  exactly  the  amount  that  the  unweighted  average  is  advanced!  The  focussing 
effect  exactly  canceled  the  Fermat  term  leaving  a  resultant  equal  to  the  geometry  effect 
alone.  This  occurs  because  a  negative  fluctuation,  which  delays  the  pulse,  acts  as  a 
converging  lens  to  increase  the  intensity. 


Ths  simple  exsmple  of  s  phase  screes  is  (he  weak-fluctuatioa  geometrical-optics 
limit  has  illustrated  our  point.  Other  remarks  on  generalisations  to  extended  media 
and  strong  fluctuations,  as  well  as  a  rigorous  extension  of  these  results  by  means  of  a 
path-integral  method  to  include  diffractive  effects  in  a  power-law  medium,  are  included 
in  Reference  1. 

There  is  no  difficulty  in  extending  the  above  results  from  a  phase  screen  to 
extended  media  in  which  (1.6)  and  (1.10)  are  replaced  by 

<  T  >  -  T*  «  -0.5C01  Jit  A~\t)\S  it*  p^t*)]  (1-11) 

<1T>  -  T0  -  +0.5CO-1  /  it  A-l(t)  [/  it*  />„(*, t*)\  (1.12) 

0  =  <dBv{t)dav(t')>  (1.13) 

where  dM  n(t)  is  the  transverse  gradient  of  the  refractive  index  due  to  the  fluctuations 
at  location  t  along  the  unperturbed  ray.  These  results  require  the  Markov  approxima¬ 
tion  (that  is,  the  quantity  in  square  brackets  in  (1.11)  is  a  local  function  of  z).  If  an 
incident  plane  wave  rather  than  a  point  source  is  used,  all  three  terms  (geometry,  Fer¬ 
mat,  and  focussing)  are  reduced  by  a  factor  of  three.  If  the  Markov  approximation  is 
not  made,  the  ratio  between  Fermat  and  geometry  remains  —2,  while  all  terms  are 

modified  by  terms  of  order  Lf/R ,  where  Lf  is  the  longitudinal  correlation  length  of  the 

medium  fluctuations. 

An  important  modification  of  the  above  result  occurs  if,  in  the  absence  of  fluctua¬ 
tions,  the  medium  has  focussing  properties.  In  ocean  sound  propagation  this  is  due  to 
the  sound  channel.  In  radio  wave  propagation  from  pulsars  this  might  be  due  to  very 
large-scale  medium  fluctuations  that  are  effectively  frozen  during  the  time  of  observa¬ 
tion.  The  modification  can  be  simply  expressed  by  generalizing  A_1(i)  in  (1.11-1.12)  for 
a  curved  ray.&  The  key  result  is  that  A~1(z)  can  be  negative  for  various  regions  along 
the  ray,  and  hence  the  geometry  term  can  be  negative  for  curvei  rays.  This  complication 
is  crucial  to  the  comparison  between  calculation  and  experiment  in  the  ocean,  though 
probably  not  in  other  media.  Note  that  this  effect  provides  another,  different  mechan¬ 
ism  by  which  fluctuations  in  a  medium  may  cause  an  average  pulse  advance. 

A  numerical  study  of  (1.11)  for  acoustic  rays  through  ocean  internal  waves  has 
revealed  this  effect  as  a  significant  bias  in  the  measurement  of  large-scale  ocean  eddies 
by  acoustic  means.  Typical  examples  of  experiments  would  be  waves  travelling  over  a 


range  of  1,000  km,  whose  travel  times  are  modified  by  about  20  ms  due  to  a  single  eddy 
of  100-km  eixe.  Tbe  effect  of  random  internal  wares  along  this  path  is  calculated  to  be 
of  tbe  same  order.  Therefore  an  average  change  in  tbe  travel  time  of  a  long-range 
ocean  ray  could  be  interpreted  either  as  the  effect  of  an  eddy  wandering  across  tbe  ray, 
or  as  the  effect  of  the  entire  internal-wave  field  changing  its  strength. 

IL  Comparison  between  Moment  Equations  and  Path-Integral  Expressions  for 
Wave  Propagation  in  Random  Media 

Many  problems  in  wave  propagation  through  random  media  concern  phenomena  in 
which  there  is  no  significant  backscatter,  so  that  a  parabolic  approximation  may  be 
made  to  the  wave  equation.6  In  these  cases  a  further  approximation,  called  the  Markov 
approximation,7  leads  to  relatively  tractable  mathematical  expressions  for  moments  of 
the  field  that  can  be  used  for  practical  calculations.  Two  quite  different  formalisms 
have  been  used  in  this  context:  the  moment-equation  and  path-integral  techniques. 

A  path-integral  expression  for  a  general  moment  of  the  field  of  a  wave  propagating 
through  an  inhomogeneous,  anisotropic  medium  in  the  presence  of  a  deterministic  back¬ 
ground  refractive  index  has  been  derived,6  and  the  expression  has  been  used  for  specific 
calculations.4, 4,9 

Moment  equations  in  coordinate  representation  have  been  derived  for  homogene¬ 
ous  isotropic  media  in  the  absence  of  a  deterministic  background.7  Treatments  of  inho¬ 
mogeneity,  anisotropy,  and  deterministic  background  by  moment-equation  techniques 
have  heretofore  been  confined  to  special  cases  involving  the  first  and  second 
moments.10,11 

We  present  here  general  moment  equations  in  coordinate  representation  that 
account  for  inhomogeneity,  anisotropy,  and  deterministic  background,  but  require  the 
Markov  approximation.  We  have  derived  these  equations1*  using  the  time-ordered- 
product  method  of  Van  Kampen,1*  which  also  provides  a  derivation  of  equations  that 
are  valid  under  conditions  more  general  than  the  Markov  approximation.  The  modified 
equations  are  more  complicated  than  those  that  require  the  Markov  approximation:  a 
special  case  was  previously  derived  by  Besieris  and  Tappert.14 
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Ws  have  also  shown  that  our  new  genera)  moment  equations  derived  under  the 
Markov  approximation  are  mathematically  equivalent  to  the  path-integral  expressions 
for  the  moments  that  have  been  previously  presented.  Thus,  the  two  popular  formal¬ 
isms,  under  the  Markov  approximation,  are  not  different  in  content. 

Consider  waves  travelling  predominantly  in  the  x  direction.  Let  f*  be  a  transverse 
coordinate  (e.g.  two-dimensional,  but  in  fact  general),  and  k  be  a  reference  wave 
number  ( k  —  2ir w/C0),  where  w  is  the  wave  frequency  and  C0  is  a  reference  wave 
speed).  Express  the  full  wave  field  as 

exp  [  £fc(*  -  C0I)]  (2.1) 

Let  the  wave  speed  (a  function  of  position  only)  be 

C(r,  x)  =  c0  [  1  -  2  t/0(r)  -  2 x)  ] "  *  «  c0  [  1  +  c/0(r)  +  rfr,  z)  ]  (2.2) 

where  U0  represents  the  deterministic  background  and  p  represents  the  fluctuating  ran¬ 
dom  medium,  assumed  to  be  a  realization  of  a  zero-mean  Gaussian  process. 

Then,  the  parabolic  equation  (in  rectangular  coordinates)  for  the  reduced  wave 
function  0  is: 

~  vty  +  *2c/o(r)0  +  *V(r,  *)0  (2.3) 

lx 

where  V*  is  the  transverse  Laplacian. 

A  moment  T  is  the  ensemble  expectation  value  of  a  product  of  0’s  and  0*’s  where 

each  0  or  0  is  evaluated  at  a  different  position  Jj  and  wavenumber  Ay.  We  write,  in 
abbreviated  form, 


r„  =  <  A*  - 

V’m  ^*+1  *”  ^m  +  n  ^ 

Define  an  operator  L0  such  that 

m+  a  j 

i.  -  E  ±-r 
1  *> 

( +*/"«) 

(25) 

The  terms  that  apply  to  the  0* s  use  the  plus  sign  and  those  that  apply  to  the  0*’s  use 
the  minus  sign.  The  subscript  j  requires  that  V*  operate  only  on  f?  and  U0j  =  U0(SJ). 


Define  the  important  combination  of  fluctuation  quantities  as 

m<«)  -  s' ±  *,*<«;.«)  (uj 

Our  general  moment  equation  under  the  Markov  approximation  can  be  written 

»,r„W  -  -  .  i0 r_M  - 1  /  it'  <  >  r..(j)  (2.7) 

where  Mtk^(zr)  is  obtained  by  evaluating  M[z)  with  all  the  Sj  at  z  shifted  by  the 
transverse  distance  that  a  deterministic  ray  through  (Sj ,  z)  moves  in  travelling  from  z 
to  z 1  (see  Figure  1).  In  other  words  Mtkyi(zr)  is  evaluated  at  point  B:  i.e.  SJ  =  x^tt(zt) 
where  the  ray  is  forced  to  go  through  (z).  The  particular  ray  is  determined  not  only 
by  the  local  position  (ij ,  z),  but  also  by  the  initial  conditions  on  the  moment;  for  exam¬ 
ple,  the  location  of  a  point  source,  or  the  direction  of  a  plane  wave.15  The  unphysical 
assumption  of  delta-correlated  medium  fluctuations  along  the  propagation  direction 
would  imply  that  {z1)  would  be  evaluated  at  point  C:  i.e.  Xj(z)  (and  z1).  In  the 
isotropic  case  (or  in  the  case  of  propagation  along  a  principal  axis  of  the  anisotropy) 
the  difference  between  evaluating  M^^z1)  at  Sj[z)  and  (*0  is  negligible,  and  the 
delta-correlated  assumption  is  adequate.  In  the  snisotropic  case,  the  necessity  of 
defining  the  unperturbed  ray  makes  (2.7)  somewhat  complicated  to  apply  for  general 
initial  conditions.  However,  since  (2.7)  is  a  linear  equation,  superposition  can  be  used 
whether  the  source  is  a  point,  an  incident  plane  wave,  or  an  arbitrary  coherent  or 
incoherent  sum  of  point  sources. 

We  now  turn  to  the  path  integral  method.  Equation  (2.3)  has  the  formal  solution 

^  (2.8) 

where  /  Di*(z)  means  integration  over  paths,  x*(z)  is  a  transverse  vector  indicating  the 
position  of  the  path  at  z ,  and 

s  “  *  1  *  ihr)  -  *>  (*•») 

In  order  to  obtain  a  given  moment,  expressions  like  (2.8)  (or  its  complex  conjugate)  are 
multiplied  together,  and  the  ensemble  average  is  taken: 


(2.10) 
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The  Markov  approximation  yields  (See  Section  V): 
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(2.11) 

We  have  shown  that  the  moment  equations  (2.7)  and  the  path  integral  expressions 
(2.11)  are  mathematically  equivalent.12 

Moment  equations  can  be  formulated  in  a  variety  of  coordinate  systems,  while 
path  integrals  require  a  rectangular  coordinate  system.  There  has  been  a  fair  amount 
of  effort  expended  on  using  polar  coordinate  systems,  especially  for  point  source  prob¬ 
lems. 

The  same  results  (for  point  sources  among  others)  can  be  obtained  in  either  polar 
or  rectangular  coordinates.  Thus,  the  results  of  Shishov12  on  the  intensity  correlation, 
derived  in  spherical  polar  coordinates,  can  be  seen  to  be  identical  (after  an  appropriate 
transformation)  to  the  results  of  Codona  et  of.,17  derived  in  rectangular  coordinates.  It 
was  necessary  for  Shishov  to  make  small  angle  approximations  in  addition  to  the  para¬ 
bolic  approximation  of  dropping  the  second  derivative  in  the  propagation  direction, 
whereas  Codona  et  at.  only  require  the  single  parabolic  approximation. 

We  have  derived  moment  equations  in  coordinate  representation  under  the  Mar¬ 
kov  approximation  that  apply  in  anisotropic,  inhomogeneous  media  with  deterministic 
background.  The  derivation  shows  the  relationship  between  these  moment  equations 
and  modified  equations  that  are  valid  under  approximations  weaker  than  Markov;  the 
second-moment  equation  of  Besieris  and  Tappert  is  a  special  case  of  these  modified 
equations. 

In  a  hierarchy  of  approximations  we  begin  with  the  parabolic  wave  equation  itself. 
A  path  integral  with  non-local  exponent  can  be  written  as  an  exact  solution,  although  it 
is  not  yet  useful  in  practice.  The  next  level  is  the  approximation  that  the  interaction 
strength  over  a  correlation  length  is  small— this  “first-order  perturbation  theory”  leads 
to  the  modified  moment  equations,  and  in  homogeneous,  isotropic  media,  to  the 


standard  moment  equations  and  path-integral  expressions.  In  anisotropic,  inhomogens* 
ous  media,  however,  a  further  approximation  is  necessary  to  obtain  the  moment  equa- 
tions  and  path  integral  expressions.  This  further  approximation  is  that  the  significant 
flow  of  wave  energy,  or  the  important  paths,  are  parallel  to  the  unperturbed  ray;  we 
call  this  the  Markov  approximation  because  its  violation  implies  the  appearance  of 
correlations  between  successive  scatterings.  We  have  shown  that  the  moment  equations 
and  the  path-integral  expressions  for  the  moments  are  mathematically  equivalent  under 
the  Markov  approximation.  Thus  the  two  formalisms  have  exactly  the  same  physical 
content.  In  an  anisotropic  medium,  the  moment  equation  involves  a  shift  operation  to 
calculate  the  medium  correlation  function  along  the  unperturbed  ray;  this  form  of  the 
moment  equation  has  not  been  given  before. 

ID.  Series  Expansion  of  the  Fourth  Moment 

We  have  developed  a  series  expression  for  the  fourth  moment  of  a  beamed  field 
incident  on  a  random  phase  screen  or  an  extended  medium.17  The  series  has  a  symmetry 
that  allows  its  first  few  terms  to  generate  useful  approximations  at  both  low  and  high 
spatial  frequency.  The  parabolic  wave  equation,  the  Markov  approximation,  and  Gaus¬ 
sian  refractive  index  fluctuations  are  assumed.  The  result  for  the  phase  screen  is 
obtained  by  Green’s-function  techniques.  The  extended-medium  result  b  derived  in  an 
analogous  manner  using  path  integral  methods.  The  same  results  can  also  be  derived 
by  moment-equation  methods.  The  behavior  of  certain  leading  terms  agree  with  previ¬ 
ous  results  for  plane-wave  and  point-source  geometries. 

We  consider  waves  propagating  from  an  arbitrary  source  distribution  in  a  random 
medium.  We  assume  the  statbtics  of  the  medium  are  locally  homogeneous,  and  we 
make  the  Markov  approximation;  i.e.  the  field  fluctuations  induced  within  a  correlation 
length  along  the  propagation  direction  are  weak.  For  a  more  complete  discussion  see 
another  of  our  papers.1*  The  wave  propagation  b  characterised  by  narrow  angular 
scattering  due  to  the  small  random  fluctuations  in  refractive  index.  It  b  then  con¬ 
venient  to  write  the  complex  monochromatic  scalar  field  as  F(Stx)ett(  where  *  b  the 
propagation  direction,  S*ia  the  transverse  coordinate  and  k  b  the  wavenumber  of  the 


wave  with  no  refractive  index  fluctuations. 

The  random  nature  of  the  fields  is  conveniently  described  by  statistical  moments 
evaluated  in  the  transverse  plane  located  at  distance  R.  Ensemble  averages  of  random 
variables  are  denoted  by  <  >.  The  first  moment 

17**)  -  <E(*,R)>  (3.1) 

or  average  of  the  field  and  the  second  moment 

r A?vr»R)  -  <E{rl,R)E*(ri,R)>  (3.2) 

or  mutual  coherence  function  are  well  understood.7  However,  there  are  few  analytic 
results  for  the  fourth  moment 

IV*;, x\,R)  -  <E(x\,R )E'(r2,R)E(?3> R)E'(x'i,R)>  (3.3) 

Previous  theoretical  work  concentrated  on  plane-wave  and  point-source  geometry.  We 
present  three  main  results  for  arbitrary  source  distribution. 

A  series  expression  for  the  fourth  moment  is  derived  as  an  expansion  of  the 
Green’s  function  for  the  fourth  moment,  thus  avoiding  the  difficulties  associated  with 
the  source  distribution.  For  the  thin-screen  problem,  the  expansion  quantity  is  a  combi¬ 
nation  of  phase  structure  functions.  For  the  extended  random  media,  the  expansion 
quantity  is  an  analogous  combination  of  phase  structure  function  densities.  The 
Green’s  function  is  expressed  as  a  multiple  path  integral.  The  resulting  series  of  path 
integrals  is  evaluated  with  a  useful  identity. 

Our  second  result  is  the  generation  of  two  series  for  the  intensity  correlation  or 
intensity  spectrum.  The  fourth  moment  r4(£J,£^f^  £},/?)  has  the  obvious  symmetries 
that  it  is  unchanged  by  interchanging  £j  and  ££  or  by  interchanging  ££  and  fj.  Each 
term  erf  the  series  expansion  does  not  share  the  symmetry  of  the  entire  expression. 
Thus  two  separate  series  are  obtained  by  invoking  symmetry.  In  principle,  either  series 
could  be  summed  to  give  T4.  We  demonstrate,  however,  that  it  is  better  to  consider 
both  series  in  order  to  describe  the  fourth  moment  with  the  fewest  number  of  terms. 
This  assertion  is  demonstrated  for  the  second  moment  of  intensity  or  intensity  correla¬ 
tion,  <7(£|, ££,#),  which  is  a  special  case  of  the  fourth  moment,  i.e. 

c{rvrt,R)  -  </(*i,/?)/(r2l/?)>  -  r<(r„s;r2>r2)  -  r<(r„r2,r2(rj)  (3.4) 

Note  that  the  symmetry  of  the  fourth  moment  has  been  explicitly  indicated.  A  clear 


presentation  of  the  behavior  of  the  intensity  correlation  series  obtained  from  the  fourth 
moment  expansion  requires  the  introduction  of  a  spatial  spectrum  of  intensity  fluctua¬ 
tions  for  a  spatially  nonstationary  random  process.  We  adopt  the  definition 

«.m*)  -  ttj  /  cwr.R)  f‘*-r *r  (3.*) 

I**)  —  op 

where 

r  =  ^-(«i+ra)  n-**  (3-6) 

(Note  the  free  format  of  the  argument  list  of  functions).  The  spectrum  has  the  pro¬ 
perty 

00 

J  *or,%R)  </?  =  C(/r,0,R)  =  </(#R)*>  (3.7) 

—  OP 

It  should  be  noted  that  the  spatial  spectrum  may  depend  on  the  centroid  ft. 

Since  there  are  two  series  for  the  intensity  correlation  there  are  also  two  series  for 
the  intensity  spectrum.  The  leading  terms  of  one  series  for  Q(jt,%R)  describe  the  small 
cf  behavior  while  the  other  series  is  valid  at  high  q!  The  rate  of  convergence  of  each 
series  provides  a  criterion  for  merging  the  two  results  to  produce  a  complete  expression 
for  the  intensity  spectrum.  In  general,  an  analogous  treatment  of  the  intensity  correla¬ 
tion  series  is  not  possible  since  the  leading  terms  of  both  series  do  not  converge  to  the 
variance  as  the  spatial  separation  approaches  zero. 

Our  third  result  is  the  demonstration  of  the  equivalence  of  path  integral  and 
moment-equation  methods.  Early  theoretical  work  on  WPRM  concentrated  on 
geometrical  optics  and  the  method  of  small  perturbations.7,1*  These  two  approaches 
were  limited  to  weak  scattering  conditions.  This  restriction  was  removed  with  the 
introduction  of  differential  equations  for  the  moments  of  the  field.19  Functional  tech¬ 
niques  of  high  energy  physics  (path  integrals  and  operator  methods)  provided  another 
point  of  view  to  WPRM.20,8  The  moment  equation  method  and  functional  techniques 
are  equivalent13  and  must  generate  identical  results  when  expansions  are  performed  in 
the  same  quantity.  This  equivalence  can  be  demonstrated  by  deriving  the  same  fourth 
moment  series  expression  using  moment-equation  methods.17 


An  example  of  (he  calculation  of  the  first  two  terms  in  the  expansion  for  a  particu¬ 
lar  case  is  shown  in  Figure  2. 

IV.Suinmary 

Descriptions  have  been  given  of  three  recent  results  in  the  theory  of  wave  propaga¬ 
tion  through  random  media.  The  first  result  is  that  pulse  arrivals  can  be  advanced  in 
time  by  the  imposition  of  a  zero-mean  random  wave-speed  fluctuation  on  the  medium. 
There  are  many  subtleties  to  the  interpretation  of  this  effect;  these  subtleties  involve 
intensity  weighting,  curved  rays,  and  strength  of  fluctuations.  The  second  result  is  that 
moment  equations  can  be  written  that  are  mathematically  equivalent  to  the  path 
integral  expressions  for  a  general  moment,  but  the  moment  equations  so  obtained  under 
the  Markov  approximation  are  not  quite  the  standard  ones  used.  The  equivalent 
moment  equations  require  knowledge  of  the  unperturbed  rays  between  source  and  the 
point  of  interest.  The  third  result  is  a  series  expansion  of  the  general  fourth  moment 
with  arbitrary  source  distribution:  a  series  expansion  whose  first  few  terms  adequately 
approximate  the  high-spatial-frequency  part  of  the  intensity  spectrum.  Previous  work 
in  the  high-spatial-frequency  regime  had  required  the  summation  of  a  large  number  of 
terms. 
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Figure  1*.  Moment-equation  expression  of  the  Markov  approximation.  The  corre¬ 
lation  should  be  taken  between  a  point  at  x  (point  P)  and  an  arbitrary 
point  at  *'  (point  A).  Instead  it  is  taken  with  the  point  B,  obtained  by 
extrapolating  along  the  unperturbed  ray  from  P.  The  assumption  of 
delta-correlated  medium  fluctuations  leads  to  the  incorrect  formulation 
of  correlations  between  points  P  and  C.  The  dashed  lines  indicate  the 
idea  of  a  scattering  as  a  function  of  angle  from  point  P. 


Figure  lb.  Path-integral  expression  of  the  Markov  approximation.  The  general 
path  at  x*  (point  A)  is  approximated  by  the  path  at  x  extrapolated 
along  the  unperturbed  ray  (point  B). 


Figure  2.  The  leading  terms  of  the  intensity  spectrum  versus  normalized  spatial 
frequency,  qRf,  where  Rj  =  ( R/k)%  is  the  Fresnel  scale.  The  curves 
are  calculated  from  expressions  given  by  Gochelashvily  and  Shishov 
[1975]  for  the  case  of  plane  waves  incident  on  a  random  phase  screen 
with  a  Kolmogorov  spectrum  of  phase  fluctuations  and  00. 

The  (— )  sign  indicates  that  $*^(?)  is  negative  at  high  frequency. 


Figure 
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